dtgsy2.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 981 行 · 第 1/5 页
HTML
981 行
Z( 3, 2 ) = ZERO
Z( 4, 2 ) = D( IS, IS )
<span class="comment">*</span><span class="comment">
</span> Z( 1, 3 ) = -B( JS, JS )
Z( 2, 3 ) = -B( JS, JSP1 )
Z( 3, 3 ) = -E( JS, JS )
Z( 4, 3 ) = -E( JS, JSP1 )
<span class="comment">*</span><span class="comment">
</span> Z( 1, 4 ) = -B( JSP1, JS )
Z( 2, 4 ) = -B( JSP1, JSP1 )
Z( 3, 4 ) = ZERO
Z( 4, 4 ) = -E( JSP1, JSP1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Set up right hand side(s)
</span><span class="comment">*</span><span class="comment">
</span> RHS( 1 ) = C( IS, JS )
RHS( 2 ) = C( IS, JSP1 )
RHS( 3 ) = F( IS, JS )
RHS( 4 ) = F( IS, JSP1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Solve Z * x = RHS
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="DGETC2.398"></a><a href="dgetc2.f.html#DGETC2.1">DGETC2</a>( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
IF( IERR.GT.0 )
$ INFO = IERR
<span class="comment">*</span><span class="comment">
</span> IF( IJOB.EQ.0 ) THEN
CALL <a name="DGESC2.403"></a><a href="dgesc2.f.html#DGESC2.1">DGESC2</a>( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
$ SCALOC )
IF( SCALOC.NE.ONE ) THEN
DO 60 K = 1, N
CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
60 CONTINUE
SCALE = SCALE*SCALOC
END IF
ELSE
CALL <a name="DLATDF.413"></a><a href="dlatdf.f.html#DLATDF.1">DLATDF</a>( IJOB, ZDIM, Z, LDZ, RHS, RDSUM,
$ RDSCAL, IPIV, JPIV )
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Unpack solution vector(s)
</span><span class="comment">*</span><span class="comment">
</span> C( IS, JS ) = RHS( 1 )
C( IS, JSP1 ) = RHS( 2 )
F( IS, JS ) = RHS( 3 )
F( IS, JSP1 ) = RHS( 4 )
<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
</span><span class="comment">*</span><span class="comment"> equation.
</span><span class="comment">*</span><span class="comment">
</span> IF( I.GT.1 ) THEN
CALL DGER( IS-1, NB, -ONE, A( 1, IS ), 1, RHS( 1 ),
$ 1, C( 1, JS ), LDC )
CALL DGER( IS-1, NB, -ONE, D( 1, IS ), 1, RHS( 1 ),
$ 1, F( 1, JS ), LDF )
END IF
IF( J.LT.Q ) THEN
CALL DAXPY( N-JE, RHS( 3 ), B( JS, JE+1 ), LDB,
$ C( IS, JE+1 ), LDC )
CALL DAXPY( N-JE, RHS( 3 ), E( JS, JE+1 ), LDE,
$ F( IS, JE+1 ), LDF )
CALL DAXPY( N-JE, RHS( 4 ), B( JSP1, JE+1 ), LDB,
$ C( IS, JE+1 ), LDC )
CALL DAXPY( N-JE, RHS( 4 ), E( JSP1, JE+1 ), LDE,
$ F( IS, JE+1 ), LDF )
END IF
<span class="comment">*</span><span class="comment">
</span> ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.1 ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Build a 4-by-4 system Z * x = RHS
</span><span class="comment">*</span><span class="comment">
</span> Z( 1, 1 ) = A( IS, IS )
Z( 2, 1 ) = A( ISP1, IS )
Z( 3, 1 ) = D( IS, IS )
Z( 4, 1 ) = ZERO
<span class="comment">*</span><span class="comment">
</span> Z( 1, 2 ) = A( IS, ISP1 )
Z( 2, 2 ) = A( ISP1, ISP1 )
Z( 3, 2 ) = D( IS, ISP1 )
Z( 4, 2 ) = D( ISP1, ISP1 )
<span class="comment">*</span><span class="comment">
</span> Z( 1, 3 ) = -B( JS, JS )
Z( 2, 3 ) = ZERO
Z( 3, 3 ) = -E( JS, JS )
Z( 4, 3 ) = ZERO
<span class="comment">*</span><span class="comment">
</span> Z( 1, 4 ) = ZERO
Z( 2, 4 ) = -B( JS, JS )
Z( 3, 4 ) = ZERO
Z( 4, 4 ) = -E( JS, JS )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Set up right hand side(s)
</span><span class="comment">*</span><span class="comment">
</span> RHS( 1 ) = C( IS, JS )
RHS( 2 ) = C( ISP1, JS )
RHS( 3 ) = F( IS, JS )
RHS( 4 ) = F( ISP1, JS )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Solve Z * x = RHS
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="DGETC2.477"></a><a href="dgetc2.f.html#DGETC2.1">DGETC2</a>( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
IF( IERR.GT.0 )
$ INFO = IERR
IF( IJOB.EQ.0 ) THEN
CALL <a name="DGESC2.481"></a><a href="dgesc2.f.html#DGESC2.1">DGESC2</a>( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
$ SCALOC )
IF( SCALOC.NE.ONE ) THEN
DO 70 K = 1, N
CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
70 CONTINUE
SCALE = SCALE*SCALOC
END IF
ELSE
CALL <a name="DLATDF.491"></a><a href="dlatdf.f.html#DLATDF.1">DLATDF</a>( IJOB, ZDIM, Z, LDZ, RHS, RDSUM,
$ RDSCAL, IPIV, JPIV )
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Unpack solution vector(s)
</span><span class="comment">*</span><span class="comment">
</span> C( IS, JS ) = RHS( 1 )
C( ISP1, JS ) = RHS( 2 )
F( IS, JS ) = RHS( 3 )
F( ISP1, JS ) = RHS( 4 )
<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
</span><span class="comment">*</span><span class="comment"> equation.
</span><span class="comment">*</span><span class="comment">
</span> IF( I.GT.1 ) THEN
CALL DGEMV( <span class="string">'N'</span>, IS-1, MB, -ONE, A( 1, IS ), LDA,
$ RHS( 1 ), 1, ONE, C( 1, JS ), 1 )
CALL DGEMV( <span class="string">'N'</span>, IS-1, MB, -ONE, D( 1, IS ), LDD,
$ RHS( 1 ), 1, ONE, F( 1, JS ), 1 )
END IF
IF( J.LT.Q ) THEN
CALL DGER( MB, N-JE, ONE, RHS( 3 ), 1,
$ B( JS, JE+1 ), LDB, C( IS, JE+1 ), LDC )
CALL DGER( MB, N-JE, ONE, RHS( 3 ), 1,
$ E( JS, JE+1 ), LDE, F( IS, JE+1 ), LDF )
END IF
<span class="comment">*</span><span class="comment">
</span> ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.2 ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Build an 8-by-8 system Z * x = RHS
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="DLASET.522"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'F'</span>, LDZ, LDZ, ZERO, ZERO, Z, LDZ )
<span class="comment">*</span><span class="comment">
</span> Z( 1, 1 ) = A( IS, IS )
Z( 2, 1 ) = A( ISP1, IS )
Z( 5, 1 ) = D( IS, IS )
<span class="comment">*</span><span class="comment">
</span> Z( 1, 2 ) = A( IS, ISP1 )
Z( 2, 2 ) = A( ISP1, ISP1 )
Z( 5, 2 ) = D( IS, ISP1 )
Z( 6, 2 ) = D( ISP1, ISP1 )
<span class="comment">*</span><span class="comment">
</span> Z( 3, 3 ) = A( IS, IS )
Z( 4, 3 ) = A( ISP1, IS )
Z( 7, 3 ) = D( IS, IS )
<span class="comment">*</span><span class="comment">
</span> Z( 3, 4 ) = A( IS, ISP1 )
Z( 4, 4 ) = A( ISP1, ISP1 )
Z( 7, 4 ) = D( IS, ISP1 )
Z( 8, 4 ) = D( ISP1, ISP1 )
<span class="comment">*</span><span class="comment">
</span> Z( 1, 5 ) = -B( JS, JS )
Z( 3, 5 ) = -B( JS, JSP1 )
Z( 5, 5 ) = -E( JS, JS )
Z( 7, 5 ) = -E( JS, JSP1 )
<span class="comment">*</span><span class="comment">
</span> Z( 2, 6 ) = -B( JS, JS )
Z( 4, 6 ) = -B( JS, JSP1 )
Z( 6, 6 ) = -E( JS, JS )
Z( 8, 6 ) = -E( JS, JSP1 )
<span class="comment">*</span><span class="comment">
</span> Z( 1, 7 ) = -B( JSP1, JS )
Z( 3, 7 ) = -B( JSP1, JSP1 )
Z( 7, 7 ) = -E( JSP1, JSP1 )
<span class="comment">*</span><span class="comment">
</span> Z( 2, 8 ) = -B( JSP1, JS )
Z( 4, 8 ) = -B( JSP1, JSP1 )
Z( 8, 8 ) = -E( JSP1, JSP1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Set up right hand side(s)
</span><span class="comment">*</span><span class="comment">
</span> K = 1
II = MB*NB + 1
DO 80 JJ = 0, NB - 1
CALL DCOPY( MB, C( IS, JS+JJ ), 1, RHS( K ), 1 )
CALL DCOPY( MB, F( IS, JS+JJ ), 1, RHS( II ), 1 )
K = K + MB
II = II + MB
80 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Solve Z * x = RHS
</span><span class="comment">*</span><span class="comment">
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?