strsen.f.html

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

            LIWMIN = MAX( 1, NN )
         ELSE IF( <a name="LSAME.307"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'N'</span> ) ) THEN
            LWMIN = MAX( 1, N )
            LIWMIN = 1
         ELSE IF( <a name="LSAME.310"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'E'</span> ) ) THEN
            LWMIN = MAX( 1, NN )
            LIWMIN = 1
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
            INFO = -15
         ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
            INFO = -17
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.EQ.0 ) THEN
         WORK( 1 ) = LWMIN
         IWORK( 1 ) = LIWMIN
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.328"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="STRSEN.328"></a><a href="strsen.f.html#STRSEN.1">STRSEN</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( M.EQ.N .OR. M.EQ.0 ) THEN
         IF( WANTS )
     $      S = ONE
         IF( WANTSP )
     $      SEP = <a name="SLANGE.340"></a><a href="slange.f.html#SLANGE.1">SLANGE</a>( <span class="string">'1'</span>, N, N, T, LDT, WORK )
         GO TO 40
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Collect the selected blocks at the top-left corner of T.
</span><span class="comment">*</span><span class="comment">
</span>      KS = 0
      PAIR = .FALSE.
      DO 20 K = 1, N
         IF( PAIR ) THEN
            PAIR = .FALSE.
         ELSE
            SWAP = SELECT( K )
            IF( K.LT.N ) THEN
               IF( T( K+1, K ).NE.ZERO ) THEN
                  PAIR = .TRUE.
                  SWAP = SWAP .OR. SELECT( K+1 )
               END IF
            END IF
            IF( SWAP ) THEN
               KS = KS + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Swap the K-th block to position KS.
</span><span class="comment">*</span><span class="comment">
</span>               IERR = 0
               KK = K
               IF( K.NE.KS )
     $            CALL <a name="STREXC.367"></a><a href="strexc.f.html#STREXC.1">STREXC</a>( COMPQ, N, T, LDT, Q, LDQ, KK, KS, WORK,
     $                         IERR )
               IF( IERR.EQ.1 .OR. IERR.EQ.2 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Blocks too close to swap: exit.
</span><span class="comment">*</span><span class="comment">
</span>                  INFO = 1
                  IF( WANTS )
     $               S = ZERO
                  IF( WANTSP )
     $               SEP = ZERO
                  GO TO 40
               END IF
               IF( PAIR )
     $            KS = KS + 1
            END IF
         END IF
   20 CONTINUE
<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">        Solve Sylvester equation for R:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           T11*R - R*T22 = scale*T12
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SLACPY.392"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'F'</span>, N1, N2, T( 1, N1+1 ), LDT, WORK, N1 )
         CALL <a name="STRSYL.393"></a><a href="strsyl.f.html#STRSYL.1">STRSYL</a>( <span class="string">'N'</span>, <span class="string">'N'</span>, -1, N1, N2, T, LDT, T( N1+1, N1+1 ),
     $                LDT, WORK, N1, SCALE, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Estimate the reciprocal of the condition number of the cluster
</span><span class="comment">*</span><span class="comment">        of eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span>         RNORM = <a name="SLANGE.399"></a><a href="slange.f.html#SLANGE.1">SLANGE</a>( <span class="string">'F'</span>, N1, N2, WORK, N1, WORK )
         IF( RNORM.EQ.ZERO ) THEN
            S = ONE
         ELSE
            S = SCALE / ( SQRT( SCALE*SCALE / RNORM+RNORM )*
     $          SQRT( RNORM ) )
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( WANTSP ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Estimate sep(T11,T22).
</span><span class="comment">*</span><span class="comment">
</span>         EST = ZERO
         KASE = 0
   30    CONTINUE
         CALL <a name="SLACN2.415"></a><a href="slacn2.f.html#SLACN2.1">SLACN2</a>( NN, WORK( NN+1 ), WORK, IWORK, EST, KASE, ISAVE )
         IF( KASE.NE.0 ) THEN
            IF( KASE.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Solve  T11*R - R*T22 = scale*X.
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="STRSYL.421"></a><a href="strsyl.f.html#STRSYL.1">STRSYL</a>( <span class="string">'N'</span>, <span class="string">'N'</span>, -1, N1, N2, T, LDT,
     $                      T( N1+1, N1+1 ), LDT, WORK, N1, SCALE,
     $                      IERR )
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Solve  T11'*R - R*T22' = scale*X.
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="STRSYL.428"></a><a href="strsyl.f.html#STRSYL.1">STRSYL</a>( <span class="string">'T'</span>, <span class="string">'T'</span>, -1, N1, N2, T, LDT,
     $                      T( N1+1, N1+1 ), LDT, WORK, N1, SCALE,
     $                      IERR )
            END IF
            GO TO 30
         END IF
<span class="comment">*</span><span class="comment">
</span>         SEP = SCALE / EST
      END IF
<span class="comment">*</span><span class="comment">
</span>   40 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Store the output eigenvalues in WR and WI.
</span><span class="comment">*</span><span class="comment">
</span>      DO 50 K = 1, N
         WR( K ) = T( K, K )
         WI( K ) = ZERO
   50 CONTINUE
      DO 60 K = 1, N - 1
         IF( T( K+1, K ).NE.ZERO ) THEN
            WI( K ) = SQRT( ABS( T( K, K+1 ) ) )*
     $                SQRT( ABS( T( K+1, K ) ) )
            WI( K+1 ) = -WI( K )
         END IF
   60 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      WORK( 1 ) = LWMIN
      IWORK( 1 ) = LIWMIN
<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="STRSEN.459"></a><a href="strsen.f.html#STRSEN.1">STRSEN</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>