dhgeqz.f.html

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

HTML
986
字号
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LWORK   (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The dimension of the array WORK.  LWORK &gt;= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          If LWORK = -1, then a workspace query is assumed; the routine
</span><span class="comment">*</span><span class="comment">          only calculates the optimal size of the WORK array, returns
</span><span class="comment">*</span><span class="comment">          this value as the first entry of the WORK array, and no error
</span><span class="comment">*</span><span class="comment">          message related to LWORK is issued by <a name="XERBLA.187"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0: successful exit
</span><span class="comment">*</span><span class="comment">          &lt; 0: if INFO = -i, the i-th argument had an illegal value
</span><span class="comment">*</span><span class="comment">          = 1,...,N: the QZ iteration did not converge.  (H,T) is not
</span><span class="comment">*</span><span class="comment">                     in Schur form, but ALPHAR(i), ALPHAI(i), and
</span><span class="comment">*</span><span class="comment">                     BETA(i), i=INFO+1,...,N should be correct.
</span><span class="comment">*</span><span class="comment">          = N+1,...,2*N: the shift calculation failed.  (H,T) is not
</span><span class="comment">*</span><span class="comment">                     in Schur form, but ALPHAR(i), ALPHAI(i), and
</span><span class="comment">*</span><span class="comment">                     BETA(i), i=INFO-N+1,...,N should be correct.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Further Details
</span><span class="comment">*</span><span class="comment">  ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Iteration counters:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  JITER  -- counts iterations.
</span><span class="comment">*</span><span class="comment">  IITER  -- counts iterations run since ILAST was last
</span><span class="comment">*</span><span class="comment">            changed.  This is therefore reset only when a 1-by-1 or
</span><span class="comment">*</span><span class="comment">            2-by-2 block deflates off the bottom.
</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><span class="comment">*</span><span class="comment">    $                     SAFETY = 1.0E+0 )
</span>      DOUBLE PRECISION   HALF, ZERO, ONE, SAFETY
      PARAMETER          ( HALF = 0.5D+0, ZERO = 0.0D+0, ONE = 1.0D+0,
     $                   SAFETY = 1.0D+2 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            ILAZR2, ILAZRO, ILPIVT, ILQ, ILSCHR, ILZ,
     $                   LQUERY
      INTEGER            ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST,
     $                   ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER,
     $                   JR, MAXIT
      DOUBLE PRECISION   A11, A12, A1I, A1R, A21, A22, A2I, A2R, AD11,
     $                   AD11L, AD12, AD12L, AD21, AD21L, AD22, AD22L,
     $                   AD32L, AN, ANORM, ASCALE, ATOL, B11, B1A, B1I,
     $                   B1R, B22, B2A, B2I, B2R, BN, BNORM, BSCALE,
     $                   BTOL, C, C11I, C11R, C12, C21, C22I, C22R, CL,
     $                   CQ, CR, CZ, ESHIFT, S, S1, S1INV, S2, SAFMAX,
     $                   SAFMIN, SCALE, SL, SQI, SQR, SR, SZI, SZR, T1,
     $                   TAU, TEMP, TEMP2, TEMPI, TEMPR, U1, U12, U12L,
     $                   U2, ULP, VS, W11, W12, W21, W22, WABS, WI, WR,
     $                   WR2
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      DOUBLE PRECISION   V( 3 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.238"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      DOUBLE PRECISION   <a name="DLAMCH.239"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLANHS.239"></a><a href="dlanhs.f.html#DLANHS.1">DLANHS</a>, <a name="DLAPY2.239"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>, <a name="DLAPY3.239"></a><a href="dlapy3.f.html#DLAPY3.1">DLAPY3</a>
      EXTERNAL           <a name="LSAME.240"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="DLAMCH.240"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLANHS.240"></a><a href="dlanhs.f.html#DLANHS.1">DLANHS</a>, <a name="DLAPY2.240"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>, <a name="DLAPY3.240"></a><a href="dlapy3.f.html#DLAPY3.1">DLAPY3</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="DLAG2.243"></a><a href="dlag2.f.html#DLAG2.1">DLAG2</a>, <a name="DLARFG.243"></a><a href="dlarfg.f.html#DLARFG.1">DLARFG</a>, <a name="DLARTG.243"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</a>, <a name="DLASET.243"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>, <a name="DLASV2.243"></a><a href="dlasv2.f.html#DLASV2.1">DLASV2</a>, DROT,
     $                   <a name="XERBLA.244"></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, DBLE, MAX, MIN, 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 JOB, COMPQ, COMPZ
</span><span class="comment">*</span><span class="comment">
</span>      IF( <a name="LSAME.253"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'E'</span> ) ) THEN
         ILSCHR = .FALSE.
         ISCHUR = 1
      ELSE IF( <a name="LSAME.256"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'S'</span> ) ) THEN
         ILSCHR = .TRUE.
         ISCHUR = 2
      ELSE
         ISCHUR = 0
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( <a name="LSAME.263"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPQ, <span class="string">'N'</span> ) ) THEN
         ILQ = .FALSE.
         ICOMPQ = 1
      ELSE IF( <a name="LSAME.266"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPQ, <span class="string">'V'</span> ) ) THEN
         ILQ = .TRUE.
         ICOMPQ = 2
      ELSE IF( <a name="LSAME.269"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPQ, <span class="string">'I'</span> ) ) THEN
         ILQ = .TRUE.
         ICOMPQ = 3
      ELSE
         ICOMPQ = 0
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( <a name="LSAME.276"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPZ, <span class="string">'N'</span> ) ) THEN
         ILZ = .FALSE.
         ICOMPZ = 1
      ELSE IF( <a name="LSAME.279"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPZ, <span class="string">'V'</span> ) ) THEN
         ILZ = .TRUE.
         ICOMPZ = 2
      ELSE IF( <a name="LSAME.282"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPZ, <span class="string">'I'</span> ) ) THEN
         ILZ = .TRUE.
         ICOMPZ = 3
      ELSE
         ICOMPZ = 0
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Check Argument Values
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
      WORK( 1 ) = MAX( 1, N )
      LQUERY = ( LWORK.EQ.-1 )
      IF( ISCHUR.EQ.0 ) THEN
         INFO = -1
      ELSE IF( ICOMPQ.EQ.0 ) THEN
         INFO = -2
      ELSE IF( ICOMPZ.EQ.0 ) THEN
         INFO = -3
      ELSE IF( N.LT.0 ) THEN
         INFO = -4
      ELSE IF( ILO.LT.1 ) THEN
         INFO = -5
      ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
         INFO = -6
      ELSE IF( LDH.LT.N ) THEN
         INFO = -8
      ELSE IF( LDT.LT.N ) THEN
         INFO = -10
      ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN
         INFO = -15
      ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN
         INFO = -17
      ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
         INFO = -19
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.318"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DHGEQZ.318"></a><a href="dhgeqz.f.html#DHGEQZ.1">DHGEQZ</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( N.LE.0 ) THEN
         WORK( 1 ) = DBLE( 1 )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Initialize Q and Z
</span><span class="comment">*</span><span class="comment">
</span>      IF( ICOMPQ.EQ.3 )
     $   CALL <a name="DLASET.334"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'Full'</span>, N, N, ZERO, ONE, Q, LDQ )
      IF( ICOMPZ.EQ.3 )
     $   CALL <a name="DLASET.336"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'Full'</span>, N, N, ZERO, ONE, Z, LDZ )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Machine Constants
</span><span class="comment">*</span><span class="comment">
</span>      IN = IHI + 1 - ILO
      SAFMIN = <a name="DLAMCH.341"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'S'</span> )
      SAFMAX = ONE / SAFMIN
      ULP = <a name="DLAMCH.343"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'E'</span> )*<a name="DLAMCH.343"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'B'</span> )
      ANORM = <a name="DLANHS.344"></a><a href="dlanhs.f.html#DLANHS.1">DLANHS</a>( <span class="string">'F'</span>, IN, H( ILO, ILO ), LDH, WORK )
      BNORM = <a name="DLANHS.345"></a><a href="dlanhs.f.html#DLANHS.1">DLANHS</a>( <span class="string">'F'</span>, IN, T( ILO, ILO ), LDT, WORK )
      ATOL = MAX( SAFMIN, ULP*ANORM )
      BTOL = MAX( SAFMIN, ULP*BNORM )
      ASCALE = ONE / MAX( SAFMIN, ANORM )
      BSCALE = ONE / MAX( SAFMIN, BNORM )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set Eigenvalues IHI+1:N
</span><span class="comment">*</span><span class="comment">
</span>      DO 30 J = IHI + 1, N
         IF( T( J, J ).LT.ZERO ) THEN
            IF( ILSCHR ) THEN
               DO 10 JR = 1, J
                  H( JR, J ) = -H( JR, J )
                  T( JR, J ) = -T( JR, J )
   10          CONTINUE
            ELSE
               H( J, J ) = -H( J, J )
               T( J, J ) = -T( J, J )
            END IF
            IF( ILZ ) THEN
               DO 20 JR = 1, N
                  Z( JR, J ) = -Z( JR, J )
   20          CONTINUE
            END IF
         END IF
         ALPHAR( J ) = H( J, J )
         ALPHAI( J ) = ZERO
         BETA( J ) = T( J, J )
   30 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If IHI &lt; ILO, skip QZ steps
</span><span class="comment">*</span><span class="comment">
</span>      IF( IHI.LT.ILO )

⌨️ 快捷键说明

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