📄 dtgsy2.f
字号:
MB = IE - IS + 1
DO 190 J = Q, P + 2, -1
*
JS = IWORK( J )
JSP1 = JS + 1
JE = IWORK( J+1 ) - 1
NB = JE - JS + 1
ZDIM = MB*NB*2
IF( ( MB.EQ.1 ) .AND. ( NB.EQ.1 ) ) THEN
*
* Build a 2-by-2 system Z' * x = RHS
*
Z( 1, 1 ) = A( IS, IS )
Z( 2, 1 ) = -B( JS, JS )
Z( 1, 2 ) = D( IS, IS )
Z( 2, 2 ) = -E( JS, JS )
*
* Set up right hand side(s)
*
RHS( 1 ) = C( IS, JS )
RHS( 2 ) = F( IS, JS )
*
* Solve Z' * x = RHS
*
CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
IF( IERR.GT.0 )
$ INFO = IERR
*
CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
IF( SCALOC.NE.ONE ) THEN
DO 130 K = 1, N
CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
130 CONTINUE
SCALE = SCALE*SCALOC
END IF
*
* Unpack solution vector(s)
*
C( IS, JS ) = RHS( 1 )
F( IS, JS ) = RHS( 2 )
*
* Substitute R(I, J) and L(I, J) into remaining
* equation.
*
IF( J.GT.P+2 ) THEN
ALPHA = RHS( 1 )
CALL DAXPY( JS-1, ALPHA, B( 1, JS ), 1, F( IS, 1 ),
$ LDF )
ALPHA = RHS( 2 )
CALL DAXPY( JS-1, ALPHA, E( 1, JS ), 1, F( IS, 1 ),
$ LDF )
END IF
IF( I.LT.P ) THEN
ALPHA = -RHS( 1 )
CALL DAXPY( M-IE, ALPHA, A( IS, IE+1 ), LDA,
$ C( IE+1, JS ), 1 )
ALPHA = -RHS( 2 )
CALL DAXPY( M-IE, ALPHA, D( IS, IE+1 ), LDD,
$ C( IE+1, JS ), 1 )
END IF
*
ELSE IF( ( MB.EQ.1 ) .AND. ( NB.EQ.2 ) ) THEN
*
* Build a 4-by-4 system Z' * x = RHS
*
Z( 1, 1 ) = A( IS, IS )
Z( 2, 1 ) = ZERO
Z( 3, 1 ) = -B( JS, JS )
Z( 4, 1 ) = -B( JSP1, JS )
*
Z( 1, 2 ) = ZERO
Z( 2, 2 ) = A( IS, IS )
Z( 3, 2 ) = -B( JS, JSP1 )
Z( 4, 2 ) = -B( JSP1, JSP1 )
*
Z( 1, 3 ) = D( IS, IS )
Z( 2, 3 ) = ZERO
Z( 3, 3 ) = -E( JS, JS )
Z( 4, 3 ) = ZERO
*
Z( 1, 4 ) = ZERO
Z( 2, 4 ) = D( IS, IS )
Z( 3, 4 ) = -E( JS, JSP1 )
Z( 4, 4 ) = -E( JSP1, JSP1 )
*
* Set up right hand side(s)
*
RHS( 1 ) = C( IS, JS )
RHS( 2 ) = C( IS, JSP1 )
RHS( 3 ) = F( IS, JS )
RHS( 4 ) = F( IS, JSP1 )
*
* Solve Z' * x = RHS
*
CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
IF( IERR.GT.0 )
$ INFO = IERR
CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
IF( SCALOC.NE.ONE ) THEN
DO 140 K = 1, N
CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
140 CONTINUE
SCALE = SCALE*SCALOC
END IF
*
* Unpack solution vector(s)
*
C( IS, JS ) = RHS( 1 )
C( IS, JSP1 ) = RHS( 2 )
F( IS, JS ) = RHS( 3 )
F( IS, JSP1 ) = RHS( 4 )
*
* Substitute R(I, J) and L(I, J) into remaining
* equation.
*
IF( J.GT.P+2 ) THEN
CALL DAXPY( JS-1, RHS( 1 ), B( 1, JS ), 1,
$ F( IS, 1 ), LDF )
CALL DAXPY( JS-1, RHS( 2 ), B( 1, JSP1 ), 1,
$ F( IS, 1 ), LDF )
CALL DAXPY( JS-1, RHS( 3 ), E( 1, JS ), 1,
$ F( IS, 1 ), LDF )
CALL DAXPY( JS-1, RHS( 4 ), E( 1, JSP1 ), 1,
$ F( IS, 1 ), LDF )
END IF
IF( I.LT.P ) THEN
CALL DGER( M-IE, NB, -ONE, A( IS, IE+1 ), LDA,
$ RHS( 1 ), 1, C( IE+1, JS ), LDC )
CALL DGER( M-IE, NB, -ONE, D( IS, IE+1 ), LDD,
$ RHS( 3 ), 1, C( IE+1, JS ), LDC )
END IF
*
ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.1 ) ) THEN
*
* Build a 4-by-4 system Z' * x = RHS
*
Z( 1, 1 ) = A( IS, IS )
Z( 2, 1 ) = A( IS, ISP1 )
Z( 3, 1 ) = -B( JS, JS )
Z( 4, 1 ) = ZERO
*
Z( 1, 2 ) = A( ISP1, IS )
Z( 2, 2 ) = A( ISP1, ISP1 )
Z( 3, 2 ) = ZERO
Z( 4, 2 ) = -B( JS, JS )
*
Z( 1, 3 ) = D( IS, IS )
Z( 2, 3 ) = D( IS, ISP1 )
Z( 3, 3 ) = -E( JS, JS )
Z( 4, 3 ) = ZERO
*
Z( 1, 4 ) = ZERO
Z( 2, 4 ) = D( ISP1, ISP1 )
Z( 3, 4 ) = ZERO
Z( 4, 4 ) = -E( JS, JS )
*
* Set up right hand side(s)
*
RHS( 1 ) = C( IS, JS )
RHS( 2 ) = C( ISP1, JS )
RHS( 3 ) = F( IS, JS )
RHS( 4 ) = F( ISP1, JS )
*
* Solve Z' * x = RHS
*
CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
IF( IERR.GT.0 )
$ INFO = IERR
*
CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
IF( SCALOC.NE.ONE ) THEN
DO 150 K = 1, N
CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
150 CONTINUE
SCALE = SCALE*SCALOC
END IF
*
* Unpack solution vector(s)
*
C( IS, JS ) = RHS( 1 )
C( ISP1, JS ) = RHS( 2 )
F( IS, JS ) = RHS( 3 )
F( ISP1, JS ) = RHS( 4 )
*
* Substitute R(I, J) and L(I, J) into remaining
* equation.
*
IF( J.GT.P+2 ) THEN
CALL DGER( MB, JS-1, ONE, RHS( 1 ), 1, B( 1, JS ),
$ 1, F( IS, 1 ), LDF )
CALL DGER( MB, JS-1, ONE, RHS( 3 ), 1, E( 1, JS ),
$ 1, F( IS, 1 ), LDF )
END IF
IF( I.LT.P ) THEN
CALL DGEMV( 'T', MB, M-IE, -ONE, A( IS, IE+1 ),
$ LDA, RHS( 1 ), 1, ONE, C( IE+1, JS ),
$ 1 )
CALL DGEMV( 'T', MB, M-IE, -ONE, D( IS, IE+1 ),
$ LDD, RHS( 3 ), 1, ONE, C( IE+1, JS ),
$ 1 )
END IF
*
ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.2 ) ) THEN
*
* Build an 8-by-8 system Z' * x = RHS
*
CALL DCOPY( LDZ*LDZ, ZERO, 0, Z, 1 )
*
Z( 1, 1 ) = A( IS, IS )
Z( 2, 1 ) = A( IS, ISP1 )
Z( 5, 1 ) = -B( JS, JS )
Z( 7, 1 ) = -B( JSP1, JS )
*
Z( 1, 2 ) = A( ISP1, IS )
Z( 2, 2 ) = A( ISP1, ISP1 )
Z( 6, 2 ) = -B( JS, JS )
Z( 8, 2 ) = -B( JSP1, JS )
*
Z( 3, 3 ) = A( IS, IS )
Z( 4, 3 ) = A( IS, ISP1 )
Z( 5, 3 ) = -B( JS, JSP1 )
Z( 7, 3 ) = -B( JSP1, JSP1 )
*
Z( 3, 4 ) = A( ISP1, IS )
Z( 4, 4 ) = A( ISP1, ISP1 )
Z( 6, 4 ) = -B( JS, JSP1 )
Z( 8, 4 ) = -B( JSP1, JSP1 )
*
Z( 1, 5 ) = D( IS, IS )
Z( 2, 5 ) = D( IS, ISP1 )
Z( 5, 5 ) = -E( JS, JS )
*
Z( 2, 6 ) = D( ISP1, ISP1 )
Z( 6, 6 ) = -E( JS, JS )
*
Z( 3, 7 ) = D( IS, IS )
Z( 4, 7 ) = D( IS, ISP1 )
Z( 5, 7 ) = -E( JS, JSP1 )
Z( 7, 7 ) = -E( JSP1, JSP1 )
*
Z( 4, 8 ) = D( ISP1, ISP1 )
Z( 6, 8 ) = -E( JS, JSP1 )
Z( 8, 8 ) = -E( JSP1, JSP1 )
*
* Set up right hand side(s)
*
K = 1
II = MB*NB + 1
DO 160 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
160 CONTINUE
*
*
* Solve Z' * x = RHS
*
CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
IF( IERR.GT.0 )
$ INFO = IERR
*
CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
IF( SCALOC.NE.ONE ) THEN
DO 170 K = 1, N
CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
170 CONTINUE
SCALE = SCALE*SCALOC
END IF
*
* Unpack solution vector(s)
*
K = 1
II = MB*NB + 1
DO 180 JJ = 0, NB - 1
CALL DCOPY( MB, RHS( K ), 1, C( IS, JS+JJ ), 1 )
CALL DCOPY( MB, RHS( II ), 1, F( IS, JS+JJ ), 1 )
K = K + MB
II = II + MB
180 CONTINUE
*
* Substitute R(I, J) and L(I, J) into remaining
* equation.
*
IF( J.GT.P+2 ) THEN
CALL DGEMM( 'N', 'T', MB, JS-1, NB, ONE,
$ C( IS, JS ), LDC, B( 1, JS ), LDB, ONE,
$ F( IS, 1 ), LDF )
CALL DGEMM( 'N', 'T', MB, JS-1, NB, ONE,
$ F( IS, JS ), LDF, E( 1, JS ), LDE, ONE,
$ F( IS, 1 ), LDF )
END IF
IF( I.LT.P ) THEN
CALL DGEMM( 'T', 'N', M-IE, NB, MB, -ONE,
$ A( IS, IE+1 ), LDA, C( IS, JS ), LDC,
$ ONE, C( IE+1, JS ), LDC )
CALL DGEMM( 'T', 'N', M-IE, NB, MB, -ONE,
$ D( IS, IE+1 ), LDD, F( IS, JS ), LDF,
$ ONE, C( IE+1, JS ), LDC )
END IF
*
END IF
*
190 CONTINUE
200 CONTINUE
*
END IF
RETURN
*
* End of DTGSY2
*
END
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -