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