zhbtrd.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 613 行 · 第 1/3 页
HTML
613 行
</span> DO 200 K = KDN + 1, 2, -1
J1 = J1 + KDN
J2 = J2 + KDN
<span class="comment">*</span><span class="comment">
</span> IF( NR.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> generate plane rotations to annihilate nonzero
</span><span class="comment">*</span><span class="comment"> elements which have been created outside the band
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="ZLARGV.401"></a><a href="zlargv.f.html#ZLARGV.1">ZLARGV</a>( NR, AB( KD1, J1-KD1 ), INCA,
$ WORK( J1 ), KD1, D( J1 ), KD1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> apply plane rotations from one side
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Dependent on the the number of diagonals either
</span><span class="comment">*</span><span class="comment"> <a name="ZLARTV.408"></a><a href="zlartv.f.html#ZLARTV.1">ZLARTV</a> or <a name="ZROT.408"></a><a href="zrot.f.html#ZROT.1">ZROT</a> is used
</span><span class="comment">*</span><span class="comment">
</span> IF( NR.GT.2*KD-1 ) THEN
DO 130 L = 1, KD - 1
CALL <a name="ZLARTV.412"></a><a href="zlartv.f.html#ZLARTV.1">ZLARTV</a>( NR, AB( KD1-L, J1-KD1+L ), INCA,
$ AB( KD1-L+1, J1-KD1+L ), INCA,
$ D( J1 ), WORK( J1 ), KD1 )
130 CONTINUE
ELSE
JEND = J1 + KD1*( NR-1 )
DO 140 JINC = J1, JEND, KD1
CALL <a name="ZROT.419"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( KDM1, AB( KD, JINC-KD ), INCX,
$ AB( KD1, JINC-KD ), INCX,
$ D( JINC ), WORK( JINC ) )
140 CONTINUE
END IF
<span class="comment">*</span><span class="comment">
</span> END IF
<span class="comment">*</span><span class="comment">
</span> IF( K.GT.2 ) THEN
IF( K.LE.N-I+1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> generate plane rotation to annihilate a(i+k-1,i)
</span><span class="comment">*</span><span class="comment"> within the band
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="ZLARTG.433"></a><a href="zlartg.f.html#ZLARTG.1">ZLARTG</a>( AB( K-1, I ), AB( K, I ),
$ D( I+K-1 ), WORK( I+K-1 ), TEMP )
AB( K-1, I ) = TEMP
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> apply rotation from the left
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="ZROT.439"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( K-3, AB( K-2, I+1 ), LDAB-1,
$ AB( K-1, I+1 ), LDAB-1, D( I+K-1 ),
$ WORK( I+K-1 ) )
END IF
NR = NR + 1
J1 = J1 - KDN - 1
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> apply plane rotations from both sides to diagonal
</span><span class="comment">*</span><span class="comment"> blocks
</span><span class="comment">*</span><span class="comment">
</span> IF( NR.GT.0 )
$ CALL <a name="ZLAR2V.451"></a><a href="zlar2v.f.html#ZLAR2V.1">ZLAR2V</a>( NR, AB( 1, J1-1 ), AB( 1, J1 ),
$ AB( 2, J1-1 ), INCA, D( J1 ),
$ WORK( J1 ), KD1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> apply plane rotations from the right
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Dependent on the the number of diagonals either
</span><span class="comment">*</span><span class="comment"> <a name="ZLARTV.459"></a><a href="zlartv.f.html#ZLARTV.1">ZLARTV</a> or <a name="ZROT.459"></a><a href="zrot.f.html#ZROT.1">ZROT</a> is used
</span><span class="comment">*</span><span class="comment">
</span> IF( NR.GT.0 ) THEN
CALL <a name="ZLACGV.462"></a><a href="zlacgv.f.html#ZLACGV.1">ZLACGV</a>( NR, WORK( J1 ), KD1 )
IF( NR.GT.2*KD-1 ) THEN
DO 150 L = 1, KD - 1
IF( J2+L.GT.N ) THEN
NRT = NR - 1
ELSE
NRT = NR
END IF
IF( NRT.GT.0 )
$ CALL <a name="ZLARTV.471"></a><a href="zlartv.f.html#ZLARTV.1">ZLARTV</a>( NRT, AB( L+2, J1-1 ), INCA,
$ AB( L+1, J1 ), INCA, D( J1 ),
$ WORK( J1 ), KD1 )
150 CONTINUE
ELSE
J1END = J1 + KD1*( NR-2 )
IF( J1END.GE.J1 ) THEN
DO 160 J1INC = J1, J1END, KD1
CALL <a name="ZROT.479"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( KDM1, AB( 3, J1INC-1 ), 1,
$ AB( 2, J1INC ), 1, D( J1INC ),
$ WORK( J1INC ) )
160 CONTINUE
END IF
LEND = MIN( KDM1, N-J2 )
LAST = J1END + KD1
IF( LEND.GT.0 )
$ CALL <a name="ZROT.487"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( LEND, AB( 3, LAST-1 ), 1,
$ AB( 2, LAST ), 1, D( LAST ),
$ WORK( LAST ) )
END IF
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span> IF( WANTQ ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> accumulate product of plane rotations in Q
</span><span class="comment">*</span><span class="comment">
</span> IF( INITQ ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> take advantage of the fact that Q was
</span><span class="comment">*</span><span class="comment"> initially the Identity matrix
</span><span class="comment">*</span><span class="comment">
</span> IQEND = MAX( IQEND, J2 )
I2 = MAX( 0, K-3 )
IQAEND = 1 + I*KD
IF( K.EQ.2 )
$ IQAEND = IQAEND + KD
IQAEND = MIN( IQAEND, IQEND )
DO 170 J = J1, J2, KD1
IBL = I - I2 / KDM1
I2 = I2 + 1
IQB = MAX( 1, J-IBL )
NQ = 1 + IQAEND - IQB
IQAEND = MIN( IQAEND+KD, IQEND )
CALL <a name="ZROT.516"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
$ 1, D( J ), WORK( J ) )
170 CONTINUE
ELSE
<span class="comment">*</span><span class="comment">
</span> DO 180 J = J1, J2, KD1
CALL <a name="ZROT.522"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
$ D( J ), WORK( J ) )
180 CONTINUE
END IF
END IF
<span class="comment">*</span><span class="comment">
</span> IF( J2+KDN.GT.N ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> adjust J2 to keep within the bounds of the matrix
</span><span class="comment">*</span><span class="comment">
</span> NR = NR - 1
J2 = J2 - KDN - 1
END IF
<span class="comment">*</span><span class="comment">
</span> DO 190 J = J1, J2, KD1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> create nonzero element a(j+kd,j-1) outside the
</span><span class="comment">*</span><span class="comment"> band and store it in WORK
</span><span class="comment">*</span><span class="comment">
</span> WORK( J+KD ) = WORK( J )*AB( KD1, J )
AB( KD1, J ) = D( J )*AB( KD1, J )
190 CONTINUE
200 CONTINUE
210 CONTINUE
END IF
<span class="comment">*</span><span class="comment">
</span> IF( KD.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> make off-diagonal elements real and copy them to E
</span><span class="comment">*</span><span class="comment">
</span> DO 220 I = 1, N - 1
T = AB( 2, I )
ABST = ABS( T )
AB( 2, I ) = ABST
E( I ) = ABST
IF( ABST.NE.ZERO ) THEN
T = T / ABST
ELSE
T = CONE
END IF
IF( I.LT.N-1 )
$ AB( 2, I+1 ) = AB( 2, I+1 )*T
IF( WANTQ ) THEN
CALL ZSCAL( N, T, Q( 1, I+1 ), 1 )
END IF
220 CONTINUE
ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> set E to zero if original matrix was diagonal
</span><span class="comment">*</span><span class="comment">
</span> DO 230 I = 1, N - 1
E( I ) = ZERO
230 CONTINUE
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> copy diagonal elements to D
</span><span class="comment">*</span><span class="comment">
</span> DO 240 I = 1, N
D( I ) = AB( 1, I )
240 CONTINUE
END IF
<span class="comment">*</span><span class="comment">
</span> RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> End of <a name="ZHBTRD.586"></a><a href="zhbtrd.f.html#ZHBTRD.1">ZHBTRD</a>
</span><span class="comment">*</span><span class="comment">
</span> END
</pre>
</body>
</html>
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?