dhgeqz.f.html

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

HTML
986 行
字号
12345
     $                 H( ILAST, ILAST-1 ), LDH, CL, SL )
            CALL DROT( ILAST+1-IFRSTM, H( IFRSTM, ILAST-1 ), 1,
     $                 H( IFRSTM, ILAST ), 1, CR, SR )
<span class="comment">*</span><span class="comment">
</span>            IF( ILAST.LT.ILASTM )
     $         CALL DROT( ILASTM-ILAST, T( ILAST-1, ILAST+1 ), LDT,
     $                    T( ILAST, ILAST+1 ), LDH, CL, SL )
            IF( IFRSTM.LT.ILAST-1 )
     $         CALL DROT( IFIRST-IFRSTM, T( IFRSTM, ILAST-1 ), 1,
     $                    T( IFRSTM, ILAST ), 1, CR, SR )
<span class="comment">*</span><span class="comment">
</span>            IF( ILQ )
     $         CALL DROT( N, Q( 1, ILAST-1 ), 1, Q( 1, ILAST ), 1, CL,
     $                    SL )
            IF( ILZ )
     $         CALL DROT( N, Z( 1, ILAST-1 ), 1, Z( 1, ILAST ), 1, CR,
     $                    SR )
<span class="comment">*</span><span class="comment">
</span>            T( ILAST-1, ILAST-1 ) = B11
            T( ILAST-1, ILAST ) = ZERO
            T( ILAST, ILAST-1 ) = ZERO
            T( ILAST, ILAST ) = B22
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           If B22 is negative, negate column ILAST
</span><span class="comment">*</span><span class="comment">
</span>            IF( B22.LT.ZERO ) THEN
               DO 210 J = IFRSTM, ILAST
                  H( J, ILAST ) = -H( J, ILAST )
                  T( J, ILAST ) = -T( J, ILAST )
  210          CONTINUE
<span class="comment">*</span><span class="comment">
</span>               IF( ILZ ) THEN
                  DO 220 J = 1, N
                     Z( J, ILAST ) = -Z( J, ILAST )
  220             CONTINUE
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Recompute shift
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="DLAG2.816"></a><a href="dlag2.f.html#DLAG2.1">DLAG2</a>( H( ILAST-1, ILAST-1 ), LDH,
     $                  T( ILAST-1, ILAST-1 ), LDT, SAFMIN*SAFETY, S1,
     $                  TEMP, WR, TEMP2, WI )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           If standardization has perturbed the shift onto real line,
</span><span class="comment">*</span><span class="comment">           do another (real single-shift) QR step.
</span><span class="comment">*</span><span class="comment">
</span>            IF( WI.EQ.ZERO )
     $         GO TO 350
            S1INV = ONE / S1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Do EISPACK (QZVAL) computation of alpha and beta
</span><span class="comment">*</span><span class="comment">
</span>            A11 = H( ILAST-1, ILAST-1 )
            A21 = H( ILAST, ILAST-1 )
            A12 = H( ILAST-1, ILAST )
            A22 = H( ILAST, ILAST )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute complex Givens rotation on right
</span><span class="comment">*</span><span class="comment">           (Assume some element of C = (sA - wB) &gt; unfl )
</span><span class="comment">*</span><span class="comment">                            __
</span><span class="comment">*</span><span class="comment">           (sA - wB) ( CZ   -SZ )
</span><span class="comment">*</span><span class="comment">                     ( SZ    CZ )
</span><span class="comment">*</span><span class="comment">
</span>            C11R = S1*A11 - WR*B11
            C11I = -WI*B11
            C12 = S1*A12
            C21 = S1*A21
            C22R = S1*A22 - WR*B22
            C22I = -WI*B22
<span class="comment">*</span><span class="comment">
</span>            IF( ABS( C11R )+ABS( C11I )+ABS( C12 ).GT.ABS( C21 )+
     $          ABS( C22R )+ABS( C22I ) ) THEN
               T1 = <a name="DLAPY3.849"></a><a href="dlapy3.f.html#DLAPY3.1">DLAPY3</a>( C12, C11R, C11I )
               CZ = C12 / T1
               SZR = -C11R / T1
               SZI = -C11I / T1
            ELSE
               CZ = <a name="DLAPY2.854"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>( C22R, C22I )
               IF( CZ.LE.SAFMIN ) THEN
                  CZ = ZERO
                  SZR = ONE
                  SZI = ZERO
               ELSE
                  TEMPR = C22R / CZ
                  TEMPI = C22I / CZ
                  T1 = <a name="DLAPY2.862"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>( CZ, C21 )
                  CZ = CZ / T1
                  SZR = -C21*TEMPR / T1
                  SZI = C21*TEMPI / T1
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute Givens rotation on left
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           (  CQ   SQ )
</span><span class="comment">*</span><span class="comment">           (  __      )  A or B
</span><span class="comment">*</span><span class="comment">           ( -SQ   CQ )
</span><span class="comment">*</span><span class="comment">
</span>            AN = ABS( A11 ) + ABS( A12 ) + ABS( A21 ) + ABS( A22 )
            BN = ABS( B11 ) + ABS( B22 )
            WABS = ABS( WR ) + ABS( WI )
            IF( S1*AN.GT.WABS*BN ) THEN
               CQ = CZ*B11
               SQR = SZR*B22
               SQI = -SZI*B22
            ELSE
               A1R = CZ*A11 + SZR*A12
               A1I = SZI*A12
               A2R = CZ*A21 + SZR*A22
               A2I = SZI*A22
               CQ = <a name="DLAPY2.887"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>( A1R, A1I )
               IF( CQ.LE.SAFMIN ) THEN
                  CQ = ZERO
                  SQR = ONE
                  SQI = ZERO
               ELSE
                  TEMPR = A1R / CQ
                  TEMPI = A1I / CQ
                  SQR = TEMPR*A2R + TEMPI*A2I
                  SQI = TEMPI*A2R - TEMPR*A2I
               END IF
            END IF
            T1 = <a name="DLAPY3.899"></a><a href="dlapy3.f.html#DLAPY3.1">DLAPY3</a>( CQ, SQR, SQI )
            CQ = CQ / T1
            SQR = SQR / T1
            SQI = SQI / T1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute diagonal elements of QBZ
</span><span class="comment">*</span><span class="comment">
</span>            TEMPR = SQR*SZR - SQI*SZI
            TEMPI = SQR*SZI + SQI*SZR
            B1R = CQ*CZ*B11 + TEMPR*B22
            B1I = TEMPI*B22
            B1A = <a name="DLAPY2.910"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>( B1R, B1I )
            B2R = CQ*CZ*B22 + TEMPR*B11
            B2I = -TEMPI*B11
            B2A = <a name="DLAPY2.913"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>( B2R, B2I )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Normalize so beta &gt; 0, and Im( alpha1 ) &gt; 0
</span><span class="comment">*</span><span class="comment">
</span>            BETA( ILAST-1 ) = B1A
            BETA( ILAST ) = B2A
            ALPHAR( ILAST-1 ) = ( WR*B1A )*S1INV
            ALPHAI( ILAST-1 ) = ( WI*B1A )*S1INV
            ALPHAR( ILAST ) = ( WR*B2A )*S1INV
            ALPHAI( ILAST ) = -( WI*B2A )*S1INV
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Step 3: Go to next block -- exit if finished.
</span><span class="comment">*</span><span class="comment">
</span>            ILAST = IFIRST - 1
            IF( ILAST.LT.ILO )
     $         GO TO 380
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Reset counters
</span><span class="comment">*</span><span class="comment">
</span>            IITER = 0
            ESHIFT = ZERO
            IF( .NOT.ILSCHR ) THEN
               ILASTM = ILAST
               IF( IFRSTM.GT.ILAST )
     $            IFRSTM = ILO
            END IF
            GO TO 350
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Usual case: 3x3 or larger block, using Francis implicit
</span><span class="comment">*</span><span class="comment">                       double-shift
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                                    2
</span><span class="comment">*</span><span class="comment">           Eigenvalue equation is  w  - c w + d = 0,
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                                         -1 2        -1
</span><span class="comment">*</span><span class="comment">           so compute 1st column of  (A B  )  - c A B   + d
</span><span class="comment">*</span><span class="comment">           using the formula in QZIT (from EISPACK)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           We assume that the block is at least 3x3
</span><span class="comment">*</span><span class="comment">
</span>            AD11 = ( ASCALE*H( ILAST-1, ILAST-1 ) ) /
     $             ( BSCALE*T( ILAST-1, ILAST-1 ) )
            AD21 = ( ASCALE*H( ILAST, ILAST-1 ) ) /
     $             ( BSCALE*T( ILAST-1, ILAST-1 ) )
            AD12 = ( ASCALE*H( ILAST-1, ILAST ) ) /
     $             ( BSCALE*T( ILAST, ILAST ) )
            AD22 = ( ASCALE*H( ILAST, ILAST ) ) /
     $             ( BSCALE*T( ILAST, ILAST ) )
            U12 = T( ILAST-1, ILAST ) / T( ILAST, ILAST )
            AD11L = ( ASCALE*H( IFIRST, IFIRST ) ) /
     $              ( BSCALE*T( IFIRST, IFIRST ) )
            AD21L = ( ASCALE*H( IFIRST+1, IFIRST ) ) /
     $              ( BSCALE*T( IFIRST, IFIRST ) )
            AD12L = ( ASCALE*H( IFIRST, IFIRST+1 )
12345

dhgeqz.f.html - 源码说明

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

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

⌨️ 快捷键说明

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