1: SUBROUTINE DGBTRS( TRANS, N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB,
2: $ INFO )
3: *
4: * -- LAPACK routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * .. Scalar Arguments ..
10: CHARACTER TRANS
11: INTEGER INFO, KL, KU, LDAB, LDB, N, NRHS
12: * ..
13: * .. Array Arguments ..
14: INTEGER IPIV( * )
15: DOUBLE PRECISION AB( LDAB, * ), B( LDB, * )
16: * ..
17: *
18: * Purpose
19: * =======
20: *
21: * DGBTRS solves a system of linear equations
22: * A * X = B or A' * X = B
23: * with a general band matrix A using the LU factorization computed
24: * by DGBTRF.
25: *
26: * Arguments
27: * =========
28: *
29: * TRANS (input) CHARACTER*1
30: * Specifies the form of the system of equations.
31: * = 'N': A * X = B (No transpose)
32: * = 'T': A'* X = B (Transpose)
33: * = 'C': A'* X = B (Conjugate transpose = Transpose)
34: *
35: * N (input) INTEGER
36: * The order of the matrix A. N >= 0.
37: *
38: * KL (input) INTEGER
39: * The number of subdiagonals within the band of A. KL >= 0.
40: *
41: * KU (input) INTEGER
42: * The number of superdiagonals within the band of A. KU >= 0.
43: *
44: * NRHS (input) INTEGER
45: * The number of right hand sides, i.e., the number of columns
46: * of the matrix B. NRHS >= 0.
47: *
48: * AB (input) DOUBLE PRECISION array, dimension (LDAB,N)
49: * Details of the LU factorization of the band matrix A, as
50: * computed by DGBTRF. U is stored as an upper triangular band
51: * matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
52: * the multipliers used during the factorization are stored in
53: * rows KL+KU+2 to 2*KL+KU+1.
54: *
55: * LDAB (input) INTEGER
56: * The leading dimension of the array AB. LDAB >= 2*KL+KU+1.
57: *
58: * IPIV (input) INTEGER array, dimension (N)
59: * The pivot indices; for 1 <= i <= N, row i of the matrix was
60: * interchanged with row IPIV(i).
61: *
62: * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
63: * On entry, the right hand side matrix B.
64: * On exit, the solution matrix X.
65: *
66: * LDB (input) INTEGER
67: * The leading dimension of the array B. LDB >= max(1,N).
68: *
69: * INFO (output) INTEGER
70: * = 0: successful exit
71: * < 0: if INFO = -i, the i-th argument had an illegal value
72: *
73: * =====================================================================
74: *
75: * .. Parameters ..
76: DOUBLE PRECISION ONE
77: PARAMETER ( ONE = 1.0D+0 )
78: * ..
79: * .. Local Scalars ..
80: LOGICAL LNOTI, NOTRAN
81: INTEGER I, J, KD, L, LM
82: * ..
83: * .. External Functions ..
84: LOGICAL LSAME
85: EXTERNAL LSAME
86: * ..
87: * .. External Subroutines ..
88: EXTERNAL DGEMV, DGER, DSWAP, DTBSV, XERBLA
89: * ..
90: * .. Intrinsic Functions ..
91: INTRINSIC MAX, MIN
92: * ..
93: * .. Executable Statements ..
94: *
95: * Test the input parameters.
96: *
97: INFO = 0
98: NOTRAN = LSAME( TRANS, 'N' )
99: IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
100: $ LSAME( TRANS, 'C' ) ) THEN
101: INFO = -1
102: ELSE IF( N.LT.0 ) THEN
103: INFO = -2
104: ELSE IF( KL.LT.0 ) THEN
105: INFO = -3
106: ELSE IF( KU.LT.0 ) THEN
107: INFO = -4
108: ELSE IF( NRHS.LT.0 ) THEN
109: INFO = -5
110: ELSE IF( LDAB.LT.( 2*KL+KU+1 ) ) THEN
111: INFO = -7
112: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
113: INFO = -10
114: END IF
115: IF( INFO.NE.0 ) THEN
116: CALL XERBLA( 'DGBTRS', -INFO )
117: RETURN
118: END IF
119: *
120: * Quick return if possible
121: *
122: IF( N.EQ.0 .OR. NRHS.EQ.0 )
123: $ RETURN
124: *
125: KD = KU + KL + 1
126: LNOTI = KL.GT.0
127: *
128: IF( NOTRAN ) THEN
129: *
130: * Solve A*X = B.
131: *
132: * Solve L*X = B, overwriting B with X.
133: *
134: * L is represented as a product of permutations and unit lower
135: * triangular matrices L = P(1) * L(1) * ... * P(n-1) * L(n-1),
136: * where each transformation L(i) is a rank-one modification of
137: * the identity matrix.
138: *
139: IF( LNOTI ) THEN
140: DO 10 J = 1, N - 1
141: LM = MIN( KL, N-J )
142: L = IPIV( J )
143: IF( L.NE.J )
144: $ CALL DSWAP( NRHS, B( L, 1 ), LDB, B( J, 1 ), LDB )
145: CALL DGER( LM, NRHS, -ONE, AB( KD+1, J ), 1, B( J, 1 ),
146: $ LDB, B( J+1, 1 ), LDB )
147: 10 CONTINUE
148: END IF
149: *
150: DO 20 I = 1, NRHS
151: *
152: * Solve U*X = B, overwriting B with X.
153: *
154: CALL DTBSV( 'Upper', 'No transpose', 'Non-unit', N, KL+KU,
155: $ AB, LDAB, B( 1, I ), 1 )
156: 20 CONTINUE
157: *
158: ELSE
159: *
160: * Solve A'*X = B.
161: *
162: DO 30 I = 1, NRHS
163: *
164: * Solve U'*X = B, overwriting B with X.
165: *
166: CALL DTBSV( 'Upper', 'Transpose', 'Non-unit', N, KL+KU, AB,
167: $ LDAB, B( 1, I ), 1 )
168: 30 CONTINUE
169: *
170: * Solve L'*X = B, overwriting B with X.
171: *
172: IF( LNOTI ) THEN
173: DO 40 J = N - 1, 1, -1
174: LM = MIN( KL, N-J )
175: CALL DGEMV( 'Transpose', LM, NRHS, -ONE, B( J+1, 1 ),
176: $ LDB, AB( KD+1, J ), 1, ONE, B( J, 1 ), LDB )
177: L = IPIV( J )
178: IF( L.NE.J )
179: $ CALL DSWAP( NRHS, B( L, 1 ), LDB, B( J, 1 ), LDB )
180: 40 CONTINUE
181: END IF
182: END IF
183: RETURN
184: *
185: * End of DGBTRS
186: *
187: END
CVSweb interface <joel.bertrand@systella.fr>