zlaqr5.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 834 行 · 第 1/5 页
HTML
834 行
$ JROW, JTOP, K, K1, KDU, KMS, KNZ, KRCOL, KZS,
$ M, M22, MBOT, MEND, MSTART, MTOP, NBMPS, NDCOL,
$ NS, NU
LOGICAL ACCUM, BLK22, BMP22
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Functions ..
</span> DOUBLE PRECISION <a name="DLAMCH.155"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>
EXTERNAL <a name="DLAMCH.156"></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"> .. Intrinsic Functions ..
</span><span class="comment">*</span><span class="comment">
</span> INTRINSIC ABS, DBLE, DCONJG, DIMAG, MAX, MIN, MOD
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Local Arrays ..
</span> COMPLEX*16 VT( 3 )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Subroutines ..
</span> EXTERNAL <a name="DLABAD.166"></a><a href="dlabad.f.html#DLABAD.1">DLABAD</a>, ZGEMM, <a name="ZLACPY.166"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>, <a name="ZLAQR1.166"></a><a href="zlaqr1.f.html#ZLAQR1.1">ZLAQR1</a>, <a name="ZLARFG.166"></a><a href="zlarfg.f.html#ZLARFG.1">ZLARFG</a>, <a name="ZLASET.166"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>,
$ ZTRMM
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Statement Functions ..
</span> DOUBLE PRECISION CABS1
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Statement Function definitions ..
</span> CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
<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"> ==== If there are no shifts, then there is nothing to do. ====
</span><span class="comment">*</span><span class="comment">
</span> IF( NSHFTS.LT.2 )
$ RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== If the active block is empty or 1-by-1, then there
</span><span class="comment">*</span><span class="comment"> . is nothing to do. ====
</span><span class="comment">*</span><span class="comment">
</span> IF( KTOP.GE.KBOT )
$ RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== NSHFTS is supposed to be even, but if is odd,
</span><span class="comment">*</span><span class="comment"> . then simply reduce it by one. ====
</span><span class="comment">*</span><span class="comment">
</span> NS = NSHFTS - MOD( NSHFTS, 2 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Machine constants for deflation ====
</span><span class="comment">*</span><span class="comment">
</span> SAFMIN = <a name="DLAMCH.195"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'SAFE MINIMUM'</span> )
SAFMAX = RONE / SAFMIN
CALL <a name="DLABAD.197"></a><a href="dlabad.f.html#DLABAD.1">DLABAD</a>( SAFMIN, SAFMAX )
ULP = <a name="DLAMCH.198"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'PRECISION'</span> )
SMLNUM = SAFMIN*( DBLE( N ) / ULP )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Use accumulated reflections to update far-from-diagonal
</span><span class="comment">*</span><span class="comment"> . entries ? ====
</span><span class="comment">*</span><span class="comment">
</span> ACCUM = ( KACC22.EQ.1 ) .OR. ( KACC22.EQ.2 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== If so, exploit the 2-by-2 block structure? ====
</span><span class="comment">*</span><span class="comment">
</span> BLK22 = ( NS.GT.2 ) .AND. ( KACC22.EQ.2 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== clear trash ====
</span><span class="comment">*</span><span class="comment">
</span> IF( KTOP+2.LE.KBOT )
$ H( KTOP+2, KTOP ) = ZERO
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== NBMPS = number of 2-shift bulges in the chain ====
</span><span class="comment">*</span><span class="comment">
</span> NBMPS = NS / 2
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== KDU = width of slab ====
</span><span class="comment">*</span><span class="comment">
</span> KDU = 6*NBMPS - 3
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Create and chase chains of NBMPS bulges ====
</span><span class="comment">*</span><span class="comment">
</span> DO 210 INCOL = 3*( 1-NBMPS ) + KTOP - 1, KBOT - 2, 3*NBMPS - 2
NDCOL = INCOL + KDU
IF( ACCUM )
$ CALL <a name="ZLASET.228"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>( <span class="string">'ALL'</span>, KDU, KDU, ZERO, ONE, U, LDU )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Near-the-diagonal bulge chase. The following loop
</span><span class="comment">*</span><span class="comment"> . performs the near-the-diagonal part of a small bulge
</span><span class="comment">*</span><span class="comment"> . multi-shift QR sweep. Each 6*NBMPS-2 column diagonal
</span><span class="comment">*</span><span class="comment"> . chunk extends from column INCOL to column NDCOL
</span><span class="comment">*</span><span class="comment"> . (including both column INCOL and column NDCOL). The
</span><span class="comment">*</span><span class="comment"> . following loop chases a 3*NBMPS column long chain of
</span><span class="comment">*</span><span class="comment"> . NBMPS bulges 3*NBMPS-2 columns to the right. (INCOL
</span><span class="comment">*</span><span class="comment"> . may be less than KTOP and and NDCOL may be greater than
</span><span class="comment">*</span><span class="comment"> . KBOT indicating phantom columns from which to chase
</span><span class="comment">*</span><span class="comment"> . bulges before they are actually introduced or to which
</span><span class="comment">*</span><span class="comment"> . to chase bulges beyond column KBOT.) ====
</span><span class="comment">*</span><span class="comment">
</span> DO 140 KRCOL = INCOL, MIN( INCOL+3*NBMPS-3, KBOT-2 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Bulges number MTOP to MBOT are active double implicit
</span><span class="comment">*</span><span class="comment"> . shift bulges. There may or may not also be small
</span><span class="comment">*</span><span class="comment"> . 2-by-2 bulge, if there is room. The inactive bulges
</span><span class="comment">*</span><span class="comment"> . (if any) must wait until the active bulges have moved
</span><span class="comment">*</span><span class="comment"> . down the diagonal to make room. The phantom matrix
</span><span class="comment">*</span><span class="comment"> . paradigm described above helps keep track. ====
</span><span class="comment">*</span><span class="comment">
</span> MTOP = MAX( 1, ( ( KTOP-1 )-KRCOL+2 ) / 3+1 )
MBOT = MIN( NBMPS, ( KBOT-KRCOL ) / 3 )
M22 = MBOT + 1
BMP22 = ( MBOT.LT.NBMPS ) .AND. ( KRCOL+3*( M22-1 ) ).EQ.
$ ( KBOT-2 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Generate reflections to chase the chain right
</span><span class="comment">*</span><span class="comment"> . one column. (The minimum value of K is KTOP-1.) ====
</span><span class="comment">*</span><span class="comment">
</span> DO 10 M = MTOP, MBOT
K = KRCOL + 3*( M-1 )
IF( K.EQ.KTOP-1 ) THEN
CALL <a name="ZLAQR1.263"></a><a href="zlaqr1.f.html#ZLAQR1.1">ZLAQR1</a>( 3, H( KTOP, KTOP ), LDH, S( 2*M-1 ),
$ S( 2*M ), V( 1, M ) )
ALPHA = V( 1, M )
CALL <a name="ZLARFG.266"></a><a href="zlarfg.f.html#ZLARFG.1">ZLARFG</a>( 3, ALPHA, V( 2, M ), 1, V( 1, M ) )
ELSE
BETA = H( K+1, K )
V( 2, M ) = H( K+2, K )
V( 3, M ) = H( K+3, K )
CALL <a name="ZLARFG.271"></a><a href="zlarfg.f.html#ZLARFG.1">ZLARFG</a>( 3, BETA, V( 2, M ), 1, V( 1, M ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== A Bulge may collapse because of vigilant
</span><span class="comment">*</span><span class="comment"> . deflation or destructive underflow. (The
</span><span class="comment">*</span><span class="comment"> . initial bulge is always collapsed.) Use
</span><span class="comment">*</span><span class="comment"> . the two-small-subdiagonals trick to try
</span><span class="comment">*</span><span class="comment"> . to get it started again. If V(2,M).NE.0 and
</span><span class="comment">*</span><span class="comment"> . V(3,M) = H(K+3,K+1) = H(K+3,K+2) = 0, then
</span><span class="comment">*</span><span class="comment"> . this bulge is collapsing into a zero
</span><span class="comment">*</span><span class="comment"> . subdiagonal. It will be restarted next
</span><span class="comment">*</span><span class="comment"> . trip through the loop.)
</span><span class="comment">*</span><span class="comment">
</span> IF( V( 1, M ).NE.ZERO .AND.
$ ( V( 3, M ).NE.ZERO .OR. ( H( K+3,
$ K+1 ).EQ.ZERO .AND. H( K+3, K+2 ).EQ.ZERO ) ) )
$ THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Typical case: not collapsed (yet). ====
</span><span class="comment">*</span><span class="comment">
</span> H( K+1, K ) = BETA
H( K+2, K ) = ZERO
H( K+3, K ) = ZERO
ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Atypical case: collapsed. Attempt to
</span><span class="comment">*</span><span class="comment"> . reintroduce ignoring H(K+1,K). If the
</span><span class="comment">*</span><span class="comment"> . fill resulting from the new reflector
</span><span class="comment">*</span><span class="comment"> . is too large, then abandon it.
</span><span class="comment">*</span><span class="comment"> . Otherwise, use the new one. ====
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="ZLAQR1.301"></a><a href="zlaqr1.f.html#ZLAQR1.1">ZLAQR1</a>( 3, H( K+1, K+1 ), LDH, S( 2*M-1 ),
$ S( 2*M ), VT )
SCL = CABS1( VT( 1 ) ) + CABS1( VT( 2 ) ) +
$ CABS1( VT( 3 ) )
IF( SCL.NE.RZERO ) THEN
VT( 1 ) = VT( 1 ) / SCL
VT( 2 ) = VT( 2 ) / SCL
VT( 3 ) = VT( 3 ) / SCL
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== The following is the traditional and
</span><span class="comment">*</span><span class="comment"> . conservative two-small-subdiagonals
</span><span class="comment">*</span><span class="comment"> . test. ====
</span><span class="comment">*</span><span class="comment"> .
</span> IF( CABS1( H( K+1, K ) )*
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?