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