1CGBSVX(1) LAPACK driver routine (version 3.1) CGBSVX(1)
2
6 CGBSVX - the LU factorization to compute the solution to a complex sys‐
7 tem of linear equations A * X = B, A**T * X = B, or A**H * X = B,
8
10 SUBROUTINE CGBSVX( FACT, TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB,
11 IPIV, EQUED, R, C, B, LDB, X, LDX, RCOND, FERR,
12 BERR, WORK, RWORK, INFO )
13 CHARACTER EQUED, FACT, TRANS
15 INTEGER INFO, KL, KU, LDAB, LDAFB, LDB, LDX, N, NRHS
17 REAL RCOND
19 INTEGER IPIV( * )
21 REAL BERR( * ), C( * ), FERR( * ), R( * ), RWORK( * )
23 COMPLEX AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ), WORK( *
25 ), X( LDX, * )
26
28 CGBSVX uses the LU factorization to compute the solution to a complex
29 system of linear equations A * X = B, A**T * X = B, or A**H * X = B,
30 where A is a band matrix of order N with KL subdiagonals and KU super‐
31 diagonals, and X and B are N-by-NRHS matrices.
32 Error bounds on the solution and a condition estimate are also pro‐
34 vided.
35
38 The following steps are performed by this subroutine:
39 1. If FACT = 'E', real scaling factors are computed to equilibrate
41 the system:
42 TRANS = 'N': diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B
43 TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
44 TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B
45 Whether or not the system will be equilibrated depends on the
46 scaling of the matrix A, but if equilibration is used, A is
47 overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')
48 or diag(C)*B (if TRANS = 'T' or 'C').
49 2. If FACT = 'N' or 'E', the LU decomposition is used to factor the
51 matrix A (after equilibration if FACT = 'E') as
52 A = L * U,
53 where L is a product of permutation and unit lower triangular
54 matrices with KL subdiagonals, and U is upper triangular with
55 KL+KU superdiagonals.
56 3. If some U(i,i)=0, so that U is exactly singular, then the routine
58 returns with INFO = i. Otherwise, the factored form of A is used
59 to estimate the condition number of the matrix A. If the
60 reciprocal of the condition number is less than machine precision,
61 INFO = N+1 is returned as a warning, but the routine still goes on
62 to solve for X and compute error bounds as described below.
63 4. The system of equations is solved for X using the factored form
65 of A.
66 5. Iterative refinement is applied to improve the computed solution
68 matrix and calculate error bounds and backward error estimates
69 for it.
70 6. If equilibration was used, the matrix X is premultiplied by
72 diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so
73 that it solves the original system before equilibration.
74
77 FACT (input) CHARACTER*1
78 Specifies whether or not the factored form of the matrix A is
79 supplied on entry, and if not, whether the matrix A should be
80 equilibrated before it is factored. = 'F': On entry, AFB and
81 IPIV contain the factored form of A. If EQUED is not 'N', the
82 matrix A has been equilibrated with scaling factors given by R
83 and C. AB, AFB, and IPIV are not modified. = 'N': The matrix
84 A will be copied to AFB and factored.
85 = 'E': The matrix A will be equilibrated if necessary, then
86 copied to AFB and factored.
87 TRANS (input) CHARACTER*1
89 Specifies the form of the system of equations. = 'N': A * X =
90 B (No transpose)
91 = 'T': A**T * X = B (Transpose)
92 = 'C': A**H * X = B (Conjugate transpose)
93 N (input) INTEGER
95 The number of linear equations, i.e., the order of the matrix
96 A. N >= 0.
97 KL (input) INTEGER
99 The number of subdiagonals within the band of A. KL >= 0.
100 KU (input) INTEGER
102 The number of superdiagonals within the band of A. KU >= 0.
103 NRHS (input) INTEGER
105 The number of right hand sides, i.e., the number of columns of
106 the matrices B and X. NRHS >= 0.
107 AB (input/output) COMPLEX array, dimension (LDAB,N)
109 On entry, the matrix A in band storage, in rows 1 to KL+KU+1.
110 The j-th column of A is stored in the j-th column of the array
111 AB as follows: AB(KU+1+i-j,j) = A(i,j) for max(1,j-
112 KU)<=i<=min(N,j+kl)
113 If FACT = 'F' and EQUED is not 'N', then A must have been equi‐
115 librated by the scaling factors in R and/or C. AB is not modi‐
116 fied if FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on
117 exit.
118 On exit, if EQUED .ne. 'N', A is scaled as follows: EQUED =
120 'R': A := diag(R) * A
121 EQUED = 'C': A := A * diag(C)
122 EQUED = 'B': A := diag(R) * A * diag(C).
123 LDAB (input) INTEGER
125 The leading dimension of the array AB. LDAB >= KL+KU+1.
126 AFB (input or output) COMPLEX array, dimension (LDAFB,N)
128 If FACT = 'F', then AFB is an input argument and on entry con‐
129 tains details of the LU factorization of the band matrix A, as
130 computed by CGBTRF. U is stored as an upper triangular band
131 matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and the
132 multipliers used during the factorization are stored in rows
133 KL+KU+2 to 2*KL+KU+1. If EQUED .ne. 'N', then AFB is the fac‐
134 tored form of the equilibrated matrix A.
135 If FACT = 'N', then AFB is an output argument and on exit
137 returns details of the LU factorization of A.
138 If FACT = 'E', then AFB is an output argument and on exit
140 returns details of the LU factorization of the equilibrated
141 matrix A (see the description of AB for the form of the equili‐
142 brated matrix).
143 LDAFB (input) INTEGER
145 The leading dimension of the array AFB. LDAFB >= 2*KL+KU+1.
146 IPIV (input or output) INTEGER array, dimension (N)
148 If FACT = 'F', then IPIV is an input argument and on entry con‐
149 tains the pivot indices from the factorization A = L*U as com‐
150 puted by CGBTRF; row i of the matrix was interchanged with row
151 IPIV(i).
152 If FACT = 'N', then IPIV is an output argument and on exit con‐
154 tains the pivot indices from the factorization A = L*U of the
155 original matrix A.
156 If FACT = 'E', then IPIV is an output argument and on exit con‐
158 tains the pivot indices from the factorization A = L*U of the
159 equilibrated matrix A.
160 EQUED (input or output) CHARACTER*1
162 Specifies the form of equilibration that was done. = 'N': No
163 equilibration (always true if FACT = 'N').
164 = 'R': Row equilibration, i.e., A has been premultiplied by
165 diag(R). = 'C': Column equilibration, i.e., A has been post‐
166 multiplied by diag(C). = 'B': Both row and column equilibra‐
167 tion, i.e., A has been replaced by diag(R) * A * diag(C).
168 EQUED is an input argument if FACT = 'F'; otherwise, it is an
169 output argument.
170 R (input or output) REAL array, dimension (N)
172 The row scale factors for A. If EQUED = 'R' or 'B', A is mul‐
173 tiplied on the left by diag(R); if EQUED = 'N' or 'C', R is not
174 accessed. R is an input argument if FACT = 'F'; otherwise, R
175 is an output argument. If FACT = 'F' and EQUED = 'R' or 'B',
176 each element of R must be positive.
177 C (input or output) REAL array, dimension (N)
179 The column scale factors for A. If EQUED = 'C' or 'B', A is
180 multiplied on the right by diag(C); if EQUED = 'N' or 'R', C is
181 not accessed. C is an input argument if FACT = 'F'; otherwise,
182 C is an output argument. If FACT = 'F' and EQUED = 'C' or 'B',
183 each element of C must be positive.
184 B (input/output) COMPLEX array, dimension (LDB,NRHS)
186 On entry, the right hand side matrix B. On exit, if EQUED =
187 'N', B is not modified; if TRANS = 'N' and EQUED = 'R' or 'B',
188 B is overwritten by diag(R)*B; if TRANS = 'T' or 'C' and EQUED
189 = 'C' or 'B', B is overwritten by diag(C)*B.
190 LDB (input) INTEGER
192 The leading dimension of the array B. LDB >= max(1,N).
193 X (output) COMPLEX array, dimension (LDX,NRHS)
195 If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
196 the original system of equations. Note that A and B are modi‐
197 fied on exit if EQUED .ne. 'N', and the solution to the equili‐
198 brated system is inv(diag(C))*X if TRANS = 'N' and EQUED = 'C'
199 or 'B', or inv(diag(R))*X if TRANS = 'T' or 'C' and EQUED = 'R'
200 or 'B'.
201 LDX (input) INTEGER
203 The leading dimension of the array X. LDX >= max(1,N).
204 RCOND (output) REAL
206 The estimate of the reciprocal condition number of the matrix A
207 after equilibration (if done). If RCOND is less than the
208 machine precision (in particular, if RCOND = 0), the matrix is
209 singular to working precision. This condition is indicated by
210 a return code of INFO > 0.
211 FERR (output) REAL array, dimension (NRHS)
213 The estimated forward error bound for each solution vector X(j)
214 (the j-th column of the solution matrix X). If XTRUE is the
215 true solution corresponding to X(j), FERR(j) is an estimated
216 upper bound for the magnitude of the largest element in (X(j) -
217 XTRUE) divided by the magnitude of the largest element in X(j).
218 The estimate is as reliable as the estimate for RCOND, and is
219 almost always a slight overestimate of the true error.
220 BERR (output) REAL array, dimension (NRHS)
222 The componentwise relative backward error of each solution vec‐
223 tor X(j) (i.e., the smallest relative change in any element of
224 A or B that makes X(j) an exact solution).
225 WORK (workspace) COMPLEX array, dimension (2*N)
227 RWORK (workspace/output) REAL array, dimension (N)
229 On exit, RWORK(1) contains the reciprocal pivot growth factor
230 norm(A)/norm(U). The "max absolute element" norm is used. If
231 RWORK(1) is much less than 1, then the stability of the LU fac‐
232 torization of the (equilibrated) matrix A could be poor. This
233 also means that the solution X, condition estimator RCOND, and
234 forward error bound FERR could be unreliable. If factorization
235 fails with 0<INFO<=N, then RWORK(1) contains the reciprocal
236 pivot growth factor for the leading INFO columns of A.
237 INFO (output) INTEGER
239 = 0: successful exit
240 < 0: if INFO = -i, the i-th argument had an illegal value
241 > 0: if INFO = i, and i is
242 <= N: U(i,i) is exactly zero. The factorization has been com‐
243 pleted, but the factor U is exactly singular, so the solution
244 and error bounds could not be computed. RCOND = 0 is returned.
245 = N+1: U is nonsingular, but RCOND is less than machine preci‐
246 sion, meaning that the matrix is singular to working precision.
247 Nevertheless, the solution and error bounds are computed
248 because there are a number of situations where the computed
249 solution can be more accurate than the value of RCOND would
250 suggest.
251 LAPACK driver routine (version 3.N1o)vember 2006 CGBSVX(1)