ctrsen.f.html

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

HTML
384 行
字号
12
</span><span class="comment">*</span><span class="comment">  I(m) is an m by m identity matrix, and kprod denotes the Kronecker
</span><span class="comment">*</span><span class="comment">  product. We estimate sigma-min(C) by the reciprocal of an estimate of
</span><span class="comment">*</span><span class="comment">  the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C)
</span><span class="comment">*</span><span class="comment">  cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  When SEP is small, small changes in T can cause large changes in
</span><span class="comment">*</span><span class="comment">  the invariant subspace. An approximate bound on the maximum angular
</span><span class="comment">*</span><span class="comment">  error in the computed right invariant subspace is
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                      EPS * norm(T) / SEP
</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
      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            LQUERY, WANTBH, WANTQ, WANTS, WANTSP
      INTEGER            IERR, K, KASE, KS, LWMIN, N1, N2, NN
      REAL               EST, RNORM, SCALE
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      INTEGER            ISAVE( 3 )
      REAL               RWORK( 1 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.201"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      REAL               <a name="CLANGE.202"></a><a href="clange.f.html#CLANGE.1">CLANGE</a>
      EXTERNAL           <a name="LSAME.203"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="CLANGE.203"></a><a href="clange.f.html#CLANGE.1">CLANGE</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="CLACN2.206"></a><a href="clacn2.f.html#CLACN2.1">CLACN2</a>, <a name="CLACPY.206"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>, <a name="CTREXC.206"></a><a href="ctrexc.f.html#CTREXC.1">CTREXC</a>, <a name="CTRSYL.206"></a><a href="ctrsyl.f.html#CTRSYL.1">CTRSYL</a>, <a name="XERBLA.206"></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, 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.215"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'B'</span> )
      WANTS = <a name="LSAME.216"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'E'</span> ) .OR. WANTBH
      WANTSP = <a name="LSAME.217"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'V'</span> ) .OR. WANTBH
      WANTQ = <a name="LSAME.218"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPQ, <span class="string">'V'</span> )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set M to the number of selected eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span>      M = 0
      DO 10 K = 1, N
         IF( SELECT( K ) )
     $      M = M + 1
   10 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      N1 = M
      N2 = N - M
      NN = N1*N2
<span class="comment">*</span><span class="comment">
</span>      INFO = 0
      LQUERY = ( LWORK.EQ.-1 )
<span class="comment">*</span><span class="comment">
</span>      IF( WANTSP ) THEN
         LWMIN = MAX( 1, 2*NN )
      ELSE IF( <a name="LSAME.237"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'N'</span> ) ) THEN
         LWMIN = 1
      ELSE IF( <a name="LSAME.239"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'E'</span> ) ) THEN
         LWMIN = MAX( 1, NN )
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( .NOT.<a name="LSAME.243"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'N'</span> ) .AND. .NOT.WANTS .AND. .NOT.WANTSP )
     $     THEN
         INFO = -1
      ELSE IF( .NOT.<a name="LSAME.246"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPQ, <span class="string">'N'</span> ) .AND. .NOT.WANTQ ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -4
      ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
         INFO = -6
      ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
         INFO = -8
      ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
         INFO = -14
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.EQ.0 ) THEN
         WORK( 1 ) = LWMIN
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.263"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CTRSEN.263"></a><a href="ctrsen.f.html#CTRSEN.1">CTRSEN</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="CLANGE.275"></a><a href="clange.f.html#CLANGE.1">CLANGE</a>( <span class="string">'1'</span>, N, N, T, LDT, RWORK )
         GO TO 40
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Collect the selected eigenvalues at the top left corner of T.
</span><span class="comment">*</span><span class="comment">
</span>      KS = 0
      DO 20 K = 1, N
         IF( SELECT( K ) ) THEN
            KS = KS + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Swap the K-th eigenvalue to position KS.
</span><span class="comment">*</span><span class="comment">
</span>            IF( K.NE.KS )
     $         CALL <a name="CTREXC.289"></a><a href="ctrexc.f.html#CTREXC.1">CTREXC</a>( COMPQ, N, T, LDT, Q, LDQ, K, KS, IERR )
         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 the 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="CLACPY.299"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'F'</span>, N1, N2, T( 1, N1+1 ), LDT, WORK, N1 )
         CALL <a name="CTRSYL.300"></a><a href="ctrsyl.f.html#CTRSYL.1">CTRSYL</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="CLANGE.306"></a><a href="clange.f.html#CLANGE.1">CLANGE</a>( <span class="string">'F'</span>, N1, N2, WORK, N1, RWORK )
         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="CLACN2.322"></a><a href="clacn2.f.html#CLACN2.1">CLACN2</a>( NN, WORK( NN+1 ), WORK, 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="CTRSYL.328"></a><a href="ctrsyl.f.html#CTRSYL.1">CTRSYL</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="CTRSYL.335"></a><a href="ctrsyl.f.html#CTRSYL.1">CTRSYL</a>( <span class="string">'C'</span>, <span class="string">'C'</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">     Copy reordered eigenvalues to W.
</span><span class="comment">*</span><span class="comment">
</span>      DO 50 K = 1, N
         W( K ) = T( K, K )
   50 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      WORK( 1 ) = LWMIN
<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="CTRSEN.357"></a><a href="ctrsen.f.html#CTRSEN.1">CTRSEN</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
12

ctrsen.f.html - 源码说明

本页面展示了「famous linear algebra library (LAPACK) ports to windows」中的 ctrsen.f.html 源码文件，采用 HTML 编程语言编写，共 384 行代码。您可以在线阅读完整代码内容，也可以返回资源详情页下载完整源码包进行本地学习和开发。

虫虫下载站收录了大量与algebra相关的技术资源，包括源代码、技术文档、电路图等，是电子工程师和嵌入式开发者的专业学习平台。

⌨️ 快捷键说明

复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?