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