ztgsen.f.html

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

HTML
678 行
字号
1234
</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 )
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+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
      DOUBLE PRECISION   DSCALE, DSUM, RDSCAL, SAFMIN
      COMPLEX*16         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>      EXTERNAL           <a name="XERBLA.342"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>, <a name="ZLACN2.342"></a><a href="zlacn2.f.html#ZLACN2.1">ZLACN2</a>, <a name="ZLACPY.342"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>, <a name="ZLASSQ.342"></a><a href="zlassq.f.html#ZLASSQ.1">ZLASSQ</a>, ZSCAL, <a name="ZTGEXC.342"></a><a href="ztgexc.f.html#ZTGEXC.1">ZTGEXC</a>,
     $                   <a name="ZTGSYL.343"></a><a href="ztgsyl.f.html#ZTGSYL.1">ZTGSYL</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, DCMPLX, DCONJG, MAX, SQRT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      DOUBLE PRECISION   <a name="DLAMCH.349"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>
      EXTERNAL           <a name="DLAMCH.350"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>
<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.374"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="ZTGSEN.374"></a><a href="ztgsen.f.html#ZTGSEN.1">ZTGSEN</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.422"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="ZTGSEN.422"></a><a href="ztgsen.f.html#ZTGSEN.1">ZTGSEN</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="ZLASSQ.439"></a><a href="zlassq.f.html#ZLASSQ.1">ZLASSQ</a>( N, A( 1, I ), 1, DSCALE, DSUM )
               CALL <a name="ZLASSQ.440"></a><a href="zlassq.f.html#ZLASSQ.1">ZLASSQ</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="DLAMCH.450"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</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="ZTGEXC.464"></a><a href="ztgexc.f.html#ZTGEXC.1">ZTGEXC</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

ztgsen.f.html - 源码说明

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

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

⌨️ 快捷键说明

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