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 >= 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"> < 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 < ILO, skip QZ steps
</span><span class="comment">*</span><span class="comment">
</span> IF( IHI.LT.ILO )
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?