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