ctgsen.f.html

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

HTML
675 行
字号
1234
</span><span class="comment">*</span><span class="comment">      Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1,
</span><span class="comment">*</span><span class="comment">      1996.
</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>      INTEGER            IDIFJB
      PARAMETER          ( IDIFJB = 3 )
      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, SWAP, WANTD, WANTD1, WANTD2, WANTP
      INTEGER            I, IERR, IJB, K, KASE, KS, LIWMIN, LWMIN, MN2,
     $                   N1, N2
      REAL               DSCALE, DSUM, RDSCAL, SAFMIN
      COMPLEX            TEMP1, TEMP2
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      INTEGER            ISAVE( 3 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      REAL               <a name="SLAMCH.342"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a> 
      EXTERNAL           <a name="CLACN2.343"></a><a href="clacn2.f.html#CLACN2.1">CLACN2</a>, <a name="CLACPY.343"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>, <a name="CLASSQ.343"></a><a href="classq.f.html#CLASSQ.1">CLASSQ</a>, CSCAL, <a name="CTGEXC.343"></a><a href="ctgexc.f.html#CTGEXC.1">CTGEXC</a>, <a name="CTGSYL.343"></a><a href="ctgsyl.f.html#CTGSYL.1">CTGSYL</a>,
     $                   <a name="SLAMCH.344"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="XERBLA.344"></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, CMPLX, CONJG, 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>      INFO = 0
      LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
<span class="comment">*</span><span class="comment">
</span>      IF( IJOB.LT.0 .OR. IJOB.GT.5 ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -5
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -7
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -9
      ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
         INFO = -13
      ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
         INFO = -15
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.371"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CTGSEN.371"></a><a href="ctgsen.f.html#CTGSEN.1">CTGSEN</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      IERR = 0
<span class="comment">*</span><span class="comment">
</span>      WANTP = IJOB.EQ.1 .OR. IJOB.GE.4
      WANTD1 = IJOB.EQ.2 .OR. IJOB.EQ.4
      WANTD2 = IJOB.EQ.3 .OR. IJOB.EQ.5
      WANTD = WANTD1 .OR. WANTD2
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set M to the dimension of the specified pair of deflating
</span><span class="comment">*</span><span class="comment">     subspaces.
</span><span class="comment">*</span><span class="comment">
</span>      M = 0
      DO 10 K = 1, N
         ALPHA( K ) = A( K, K )
         BETA( K ) = B( K, K )
         IF( K.LT.N ) THEN
            IF( SELECT( K ) )
     $         M = M + 1
         ELSE
            IF( SELECT( N ) )
     $         M = M + 1
         END IF
   10 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      IF( IJOB.EQ.1 .OR. IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN
         LWMIN = MAX( 1, 2*M*(N-M) )
         LIWMIN = MAX( 1, N+2 )
      ELSE IF( IJOB.EQ.3 .OR. IJOB.EQ.5 ) THEN
         LWMIN = MAX( 1, 4*M*(N-M) )
         LIWMIN = MAX( 1, 2*M*(N-M), N+2 )
      ELSE
         LWMIN = 1
         LIWMIN = 1
      END IF
<span class="comment">*</span><span class="comment">
</span>      WORK( 1 ) = LWMIN
      IWORK( 1 ) = LIWMIN
<span class="comment">*</span><span class="comment">
</span>      IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
         INFO = -21
      ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
         INFO = -23
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.419"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CTGSEN.419"></a><a href="ctgsen.f.html#CTGSEN.1">CTGSEN</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( WANTP ) THEN
            PL = ONE
            PR = ONE
         END IF
         IF( WANTD ) THEN
            DSCALE = ZERO
            DSUM = ONE
            DO 20 I = 1, N
               CALL <a name="CLASSQ.436"></a><a href="classq.f.html#CLASSQ.1">CLASSQ</a>( N, A( 1, I ), 1, DSCALE, DSUM )
               CALL <a name="CLASSQ.437"></a><a href="classq.f.html#CLASSQ.1">CLASSQ</a>( N, B( 1, I ), 1, DSCALE, DSUM )
   20       CONTINUE
            DIF( 1 ) = DSCALE*SQRT( DSUM )
            DIF( 2 ) = DIF( 1 )
         END IF
         GO TO 70
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Get machine constant
</span><span class="comment">*</span><span class="comment">
</span>      SAFMIN = <a name="SLAMCH.447"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'S'</span> )
<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 (A, B).
</span><span class="comment">*</span><span class="comment">
</span>      KS = 0
      DO 30 K = 1, N
         SWAP = SELECT( K )
         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. Compute unitary Q
</span><span class="comment">*</span><span class="comment">           and Z that will swap adjacent diagonal blocks in (A, B).
</span><span class="comment">*</span><span class="comment">
</span>            IF( K.NE.KS )
     $         CALL <a name="CTGEXC.461"></a><a href="ctgexc.f.html#CTGEXC.1">CTGEXC</a>( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
     $                      LDZ, K, KS, IERR )
<span class="comment">*</span><span class="comment">
</span>            IF( IERR.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Swap is rejected: exit.
</span><span class="comment">*</span><span class="comment">
</span>               INFO = 1
               IF( WANTP ) THEN
                  PL = ZERO
                  PR = ZERO
               END IF
               IF( WANTD ) THEN
                  DIF( 1 ) = ZERO
                  DIF( 2 ) = ZERO
               END IF
               GO TO 70
            END IF
         END IF
   30 CONTINUE
      IF( WANTP ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Solve generalized Sylvester equation for R and L:
</span><span class="comment">*</span><span class="comment">                   A11 * R - L * A22 = A12
</span><span class="comment">*</span><span class="comment">                   B11 * R - L * B22 = B12
</span><span class="comment">*</span><span class="comment">
</span>         N1 = M
         N2 = N - M
1234

ctgsen.f.html - 源码说明

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

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

⌨️ 快捷键说明

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