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**T * 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**T* X = B (Transpose)
33 * = 'C': A**T* 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**T*X = B.
161 *
162 DO 30 I = 1, NRHS
163 *
164 * Solve U**T*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**T*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
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**T * 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**T* X = B (Transpose)
33 * = 'C': A**T* 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**T*X = B.
161 *
162 DO 30 I = 1, NRHS
163 *
164 * Solve U**T*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**T*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