zhbgst.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 1,016 行 · 第 1/5 页
HTML
1,016 行
$ RWORK( I-K+KA-M )*AB( KA1, I-K )
RA1 = RA
END IF
END IF
J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
NR = ( N-J2+KA ) / KA1
J1 = J2 + ( NR-1 )*KA1
IF( UPDATE ) THEN
J2T = MAX( J2, I+2*KA-K+1 )
ELSE
J2T = J2
END IF
NRT = ( N-J2T+KA ) / KA1
DO 320 J = J2T, J1, KA1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> create nonzero element a(j+1,j-ka) outside the band
</span><span class="comment">*</span><span class="comment"> and store it in WORK(j-m)
</span><span class="comment">*</span><span class="comment">
</span> WORK( J-M ) = WORK( J-M )*AB( KA1, J-KA+1 )
AB( KA1, J-KA+1 ) = RWORK( J-M )*AB( KA1, J-KA+1 )
320 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> generate rotations in 1st set to annihilate elements which
</span><span class="comment">*</span><span class="comment"> have been created outside the band
</span><span class="comment">*</span><span class="comment">
</span> IF( NRT.GT.0 )
$ CALL <a name="ZLARGV.620"></a><a href="zlargv.f.html#ZLARGV.1">ZLARGV</a>( NRT, AB( KA1, J2T-KA ), INCA, WORK( J2T-M ),
$ KA1, RWORK( J2T-M ), KA1 )
IF( NR.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> apply rotations in 1st set from the left
</span><span class="comment">*</span><span class="comment">
</span> DO 330 L = 1, KA - 1
CALL <a name="ZLARTV.627"></a><a href="zlartv.f.html#ZLARTV.1">ZLARTV</a>( NR, AB( L+1, J2-L ), INCA,
$ AB( L+2, J2-L ), INCA, RWORK( J2-M ),
$ WORK( J2-M ), KA1 )
330 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> apply rotations in 1st set from both sides to diagonal
</span><span class="comment">*</span><span class="comment"> blocks
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="ZLAR2V.635"></a><a href="zlar2v.f.html#ZLAR2V.1">ZLAR2V</a>( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
$ INCA, RWORK( J2-M ), WORK( J2-M ), KA1 )
<span class="comment">*</span><span class="comment">
</span> CALL <a name="ZLACGV.638"></a><a href="zlacgv.f.html#ZLACGV.1">ZLACGV</a>( NR, WORK( J2-M ), KA1 )
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> start applying rotations in 1st set from the right
</span><span class="comment">*</span><span class="comment">
</span> DO 340 L = KA - 1, KB - K + 1, -1
NRT = ( N-J2+L ) / KA1
IF( NRT.GT.0 )
$ CALL <a name="ZLARTV.646"></a><a href="zlartv.f.html#ZLARTV.1">ZLARTV</a>( NRT, AB( KA1-L+1, J2 ), INCA,
$ AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ),
$ WORK( J2-M ), KA1 )
340 CONTINUE
<span class="comment">*</span><span class="comment">
</span> IF( WANTX ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> post-multiply X by product of rotations in 1st set
</span><span class="comment">*</span><span class="comment">
</span> DO 350 J = J2, J1, KA1
CALL <a name="ZROT.656"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
$ RWORK( J-M ), WORK( J-M ) )
350 CONTINUE
END IF
360 CONTINUE
<span class="comment">*</span><span class="comment">
</span> IF( UPDATE ) THEN
IF( I2.LE.N .AND. KBT.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> create nonzero element a(i-kbt+ka+1,i-kbt) outside the
</span><span class="comment">*</span><span class="comment"> band and store it in WORK(i-kbt)
</span><span class="comment">*</span><span class="comment">
</span> WORK( I-KBT ) = -BB( KBT+1, I-KBT )*RA1
END IF
END IF
<span class="comment">*</span><span class="comment">
</span> DO 400 K = KB, 1, -1
IF( UPDATE ) THEN
J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1
ELSE
J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> finish applying rotations in 2nd set from the right
</span><span class="comment">*</span><span class="comment">
</span> DO 370 L = KB - K, 1, -1
NRT = ( N-J2+KA+L ) / KA1
IF( NRT.GT.0 )
$ CALL <a name="ZLARTV.684"></a><a href="zlartv.f.html#ZLARTV.1">ZLARTV</a>( NRT, AB( KA1-L+1, J2-KA ), INCA,
$ AB( KA1-L, J2-KA+1 ), INCA,
$ RWORK( J2-KA ), WORK( J2-KA ), KA1 )
370 CONTINUE
NR = ( N-J2+KA ) / KA1
J1 = J2 + ( NR-1 )*KA1
DO 380 J = J1, J2, -KA1
WORK( J ) = WORK( J-KA )
RWORK( J ) = RWORK( J-KA )
380 CONTINUE
DO 390 J = J2, J1, KA1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> create nonzero element a(j+1,j-ka) outside the band
</span><span class="comment">*</span><span class="comment"> and store it in WORK(j)
</span><span class="comment">*</span><span class="comment">
</span> WORK( J ) = WORK( J )*AB( KA1, J-KA+1 )
AB( KA1, J-KA+1 ) = RWORK( J )*AB( KA1, J-KA+1 )
390 CONTINUE
IF( UPDATE ) THEN
IF( I-K.LT.N-KA .AND. K.LE.KBT )
$ WORK( I-K+KA ) = WORK( I-K )
END IF
400 CONTINUE
<span class="comment">*</span><span class="comment">
</span> DO 440 K = KB, 1, -1
J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
NR = ( N-J2+KA ) / KA1
J1 = J2 + ( NR-1 )*KA1
IF( NR.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> generate rotations in 2nd set to annihilate elements
</span><span class="comment">*</span><span class="comment"> which have been created outside the band
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="ZLARGV.717"></a><a href="zlargv.f.html#ZLARGV.1">ZLARGV</a>( NR, AB( KA1, J2-KA ), INCA, WORK( J2 ), KA1,
$ RWORK( J2 ), KA1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> apply rotations in 2nd set from the left
</span><span class="comment">*</span><span class="comment">
</span> DO 410 L = 1, KA - 1
CALL <a name="ZLARTV.723"></a><a href="zlartv.f.html#ZLARTV.1">ZLARTV</a>( NR, AB( L+1, J2-L ), INCA,
$ AB( L+2, J2-L ), INCA, RWORK( J2 ),
$ WORK( J2 ), KA1 )
410 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> apply rotations in 2nd set from both sides to diagonal
</span><span class="comment">*</span><span class="comment"> blocks
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="ZLAR2V.731"></a><a href="zlar2v.f.html#ZLAR2V.1">ZLAR2V</a>( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
$ INCA, RWORK( J2 ), WORK( J2 ), KA1 )
<span class="comment">*</span><span class="comment">
</span> CALL <a name="ZLACGV.734"></a><a href="zlacgv.f.html#ZLACGV.1">ZLACGV</a>( NR, WORK( J2 ), KA1 )
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> start applying rotations in 2nd set from the right
</span><span class="comment">*</span><span class="comment">
</span> DO 420 L = KA - 1, KB - K + 1, -1
NRT = ( N-J2+L ) / KA1
IF( NRT.GT.0 )
$ CALL <a name="ZLARTV.742"></a><a href="zlartv.f.html#ZLARTV.1">ZLARTV</a>( NRT, AB( KA1-L+1, J2 ), INCA,
$ AB( KA1-L, J2+1 ), INCA, RWORK( J2 ),
$ WORK( J2 ), KA1 )
420 CONTINUE
<span class="comment">*</span><span class="comment">
</span> IF( WANTX ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> post-multiply X by product of rotations in 2nd set
</span><span class="comment">*</span><span class="comment">
</span> DO 430 J = J2, J1, KA1
CALL <a name="ZROT.752"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
$ RWORK( J ), WORK( J ) )
430 CONTINUE
END IF
440 CONTINUE
<span class="comment">*</span><span class="comment">
</span> DO 460 K = 1, KB - 1
J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> finish applying rotations in 1st set from the right
</span><span class="comment">*</span><span class="comment">
</span> DO 450 L = KB - K, 1, -1
NRT = ( N-J2+L ) / KA1
IF( NRT.GT.0 )
$ CALL <a name="ZLARTV.766"></a><a href="zlartv.f.html#ZLARTV.1">ZLARTV</a>( NRT, AB( KA1-L+1, J2 ), INCA,
$ AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ),
$ WORK( J2-M ), KA1 )
450 CONTINUE
460 CONTINUE
<span class="comment">*</span><span class="comment">
</span> IF( KB.GT.1 ) THEN
DO 470 J = N - 1, I2 + KA, -1
RWORK( J-M ) = RWORK( J-KA-M )
WORK( J-M ) = WORK( J-KA-M )
470 CONTINUE
END IF
<span class="comment">*</span><span class="comment">
</span> END IF
<span class="comment">*</span><span class="comment">
</span> GO TO 10
<span class="comment">*</span><span class="comment">
</span> 480 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> **************************** Phase 2 *****************************
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> The logical structure of this phase is:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> UPDATE = .TRUE.
</span><span class="comment">*</span><span class="comment"> DO I = 1, M
</span><span class="comment">*</span><span class="comment"> use S(i) to update A and create a new bulge
</span><span class="comment">*</span><span class="comment"> apply rotations to push all bulges KA positions upward
</span><span class="comment">*</span><span class="comment"> END DO
</span><span class="comment">*</span><span class="comment"> UPDATE = .FALSE.
</span><span class="comment">*</span><span class="comment"> DO I = M - KA - 1, 2, -1
</span><span class="comment">*</span><span class="comment"> apply rotations to push all bulges KA positions upward
</span><span class="comment">*</span><span class="comment"> END DO
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?