ztgsyl.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 599 行 · 第 1/3 页
HTML
599 行
</span> 70 CONTINUE
IWORK( Q+1 ) = N + 1
IF( IWORK( Q ).EQ.IWORK( Q+1 ) )
$ Q = Q - 1
<span class="comment">*</span><span class="comment">
</span> IF( NOTRAN ) THEN
DO 150 IROUND = 1, ISOLVE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Solve (I, J) - subsystem
</span><span class="comment">*</span><span class="comment"> A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
</span><span class="comment">*</span><span class="comment"> D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
</span><span class="comment">*</span><span class="comment"> for I = P, P - 1, ..., 1; J = 1, 2, ..., Q
</span><span class="comment">*</span><span class="comment">
</span> PQ = 0
SCALE = ONE
DSCALE = ZERO
DSUM = ONE
DO 130 J = P + 2, Q
JS = IWORK( J )
JE = IWORK( J+1 ) - 1
NB = JE - JS + 1
DO 120 I = P, 1, -1
IS = IWORK( I )
IE = IWORK( I+1 ) - 1
MB = IE - IS + 1
CALL <a name="ZTGSY2.407"></a><a href="ztgsy2.f.html#ZTGSY2.1">ZTGSY2</a>( TRANS, IFUNC, MB, NB, A( IS, IS ), LDA,
$ B( JS, JS ), LDB, C( IS, JS ), LDC,
$ D( IS, IS ), LDD, E( JS, JS ), LDE,
$ F( IS, JS ), LDF, SCALOC, DSUM, DSCALE,
$ LINFO )
IF( LINFO.GT.0 )
$ INFO = LINFO
PQ = PQ + MB*NB
IF( SCALOC.NE.ONE ) THEN
DO 80 K = 1, JS - 1
CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ),
$ C( 1, K ), 1 )
CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ),
$ F( 1, K ), 1 )
80 CONTINUE
DO 90 K = JS, JE
CALL ZSCAL( IS-1, DCMPLX( SCALOC, ZERO ),
$ C( 1, K ), 1 )
CALL ZSCAL( IS-1, DCMPLX( SCALOC, ZERO ),
$ F( 1, K ), 1 )
90 CONTINUE
DO 100 K = JS, JE
CALL ZSCAL( M-IE, DCMPLX( SCALOC, ZERO ),
$ C( IE+1, K ), 1 )
CALL ZSCAL( M-IE, DCMPLX( SCALOC, ZERO ),
$ F( IE+1, K ), 1 )
100 CONTINUE
DO 110 K = JE + 1, N
CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ),
$ C( 1, K ), 1 )
CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ),
$ F( 1, K ), 1 )
110 CONTINUE
SCALE = SCALE*SCALOC
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Substitute R(I,J) and L(I,J) into remaining equation.
</span><span class="comment">*</span><span class="comment">
</span> IF( I.GT.1 ) THEN
CALL ZGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, IS-1, NB, MB,
$ DCMPLX( -ONE, ZERO ), A( 1, IS ), LDA,
$ C( IS, JS ), LDC, DCMPLX( ONE, ZERO ),
$ C( 1, JS ), LDC )
CALL ZGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, IS-1, NB, MB,
$ DCMPLX( -ONE, ZERO ), D( 1, IS ), LDD,
$ C( IS, JS ), LDC, DCMPLX( ONE, ZERO ),
$ F( 1, JS ), LDF )
END IF
IF( J.LT.Q ) THEN
CALL ZGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, MB, N-JE, NB,
$ DCMPLX( ONE, ZERO ), F( IS, JS ), LDF,
$ B( JS, JE+1 ), LDB,
$ DCMPLX( ONE, ZERO ), C( IS, JE+1 ),
$ LDC )
CALL ZGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, MB, N-JE, NB,
$ DCMPLX( ONE, ZERO ), F( IS, JS ), LDF,
$ E( JS, JE+1 ), LDE,
$ DCMPLX( ONE, ZERO ), F( IS, JE+1 ),
$ LDF )
END IF
120 CONTINUE
130 CONTINUE
IF( DSCALE.NE.ZERO ) THEN
IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN
DIF = SQRT( DBLE( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) )
ELSE
DIF = SQRT( DBLE( PQ ) ) / ( DSCALE*SQRT( DSUM ) )
END IF
END IF
IF( ISOLVE.EQ.2 .AND. IROUND.EQ.1 ) THEN
IF( NOTRAN ) THEN
IFUNC = IJOB
END IF
SCALE2 = SCALE
CALL <a name="ZLACPY.481"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>( <span class="string">'F'</span>, M, N, C, LDC, WORK, M )
CALL <a name="ZLACPY.482"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>( <span class="string">'F'</span>, M, N, F, LDF, WORK( M*N+1 ), M )
CALL <a name="ZLASET.483"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>( <span class="string">'F'</span>, M, N, CZERO, CZERO, C, LDC )
CALL <a name="ZLASET.484"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>( <span class="string">'F'</span>, M, N, CZERO, CZERO, F, LDF )
ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN
CALL <a name="ZLACPY.486"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>( <span class="string">'F'</span>, M, N, WORK, M, C, LDC )
CALL <a name="ZLACPY.487"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>( <span class="string">'F'</span>, M, N, WORK( M*N+1 ), M, F, LDF )
SCALE = SCALE2
END IF
150 CONTINUE
ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Solve transposed (I, J)-subsystem
</span><span class="comment">*</span><span class="comment"> A(I, I)' * R(I, J) + D(I, I)' * L(I, J) = C(I, J)
</span><span class="comment">*</span><span class="comment"> R(I, J) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
</span><span class="comment">*</span><span class="comment"> for I = 1,2,..., P; J = Q, Q-1,..., 1
</span><span class="comment">*</span><span class="comment">
</span> SCALE = ONE
DO 210 I = 1, P
IS = IWORK( I )
IE = IWORK( I+1 ) - 1
MB = IE - IS + 1
DO 200 J = Q, P + 2, -1
JS = IWORK( J )
JE = IWORK( J+1 ) - 1
NB = JE - JS + 1
CALL <a name="ZTGSY2.507"></a><a href="ztgsy2.f.html#ZTGSY2.1">ZTGSY2</a>( TRANS, IFUNC, MB, NB, A( IS, IS ), LDA,
$ B( JS, JS ), LDB, C( IS, JS ), LDC,
$ D( IS, IS ), LDD, E( JS, JS ), LDE,
$ F( IS, JS ), LDF, SCALOC, DSUM, DSCALE,
$ LINFO )
IF( LINFO.GT.0 )
$ INFO = LINFO
IF( SCALOC.NE.ONE ) THEN
DO 160 K = 1, JS - 1
CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), C( 1, K ),
$ 1 )
CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), F( 1, K ),
$ 1 )
160 CONTINUE
DO 170 K = JS, JE
CALL ZSCAL( IS-1, DCMPLX( SCALOC, ZERO ),
$ C( 1, K ), 1 )
CALL ZSCAL( IS-1, DCMPLX( SCALOC, ZERO ),
$ F( 1, K ), 1 )
170 CONTINUE
DO 180 K = JS, JE
CALL ZSCAL( M-IE, DCMPLX( SCALOC, ZERO ),
$ C( IE+1, K ), 1 )
CALL ZSCAL( M-IE, DCMPLX( SCALOC, ZERO ),
$ F( IE+1, K ), 1 )
180 CONTINUE
DO 190 K = JE + 1, N
CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), C( 1, K ),
$ 1 )
CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), F( 1, K ),
$ 1 )
190 CONTINUE
SCALE = SCALE*SCALOC
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Substitute R(I,J) and L(I,J) into remaining equation.
</span><span class="comment">*</span><span class="comment">
</span> IF( J.GT.P+2 ) THEN
CALL ZGEMM( <span class="string">'N'</span>, <span class="string">'C'</span>, MB, JS-1, NB,
$ DCMPLX( ONE, ZERO ), C( IS, JS ), LDC,
$ B( 1, JS ), LDB, DCMPLX( ONE, ZERO ),
$ F( IS, 1 ), LDF )
CALL ZGEMM( <span class="string">'N'</span>, <span class="string">'C'</span>, MB, JS-1, NB,
$ DCMPLX( ONE, ZERO ), F( IS, JS ), LDF,
$ E( 1, JS ), LDE, DCMPLX( ONE, ZERO ),
$ F( IS, 1 ), LDF )
END IF
IF( I.LT.P ) THEN
CALL ZGEMM( <span class="string">'C'</span>, <span class="string">'N'</span>, M-IE, NB, MB,
$ DCMPLX( -ONE, ZERO ), A( IS, IE+1 ), LDA,
$ C( IS, JS ), LDC, DCMPLX( ONE, ZERO ),
$ C( IE+1, JS ), LDC )
CALL ZGEMM( <span class="string">'C'</span>, <span class="string">'N'</span>, M-IE, NB, MB,
$ DCMPLX( -ONE, ZERO ), D( IS, IE+1 ), LDD,
$ F( IS, JS ), LDF, DCMPLX( ONE, ZERO ),
$ C( IE+1, JS ), LDC )
END IF
200 CONTINUE
210 CONTINUE
END IF
<span class="comment">*</span><span class="comment">
</span> WORK( 1 ) = LWMIN
<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="ZTGSYL.572"></a><a href="ztgsyl.f.html#ZTGSYL.1">ZTGSYL</a>
</span><span class="comment">*</span><span class="comment">
</span> END
</pre>
</body>
</html>
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?