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