1DPBSVX(1) LAPACK driver routine (version 3.1) DPBSVX(1)
2
6 DPBSVX - the Cholesky factorization A = U**T*U or A = L*L**T to compute
7 the solution to a real system of linear equations A * X = B,
8
10 SUBROUTINE DPBSVX( FACT, UPLO, N, KD, NRHS, AB, LDAB, AFB, LDAFB,
11 EQUED, S, B, LDB, X, LDX, RCOND, FERR, BERR, WORK,
12 IWORK, INFO )
13 CHARACTER EQUED, FACT, UPLO
15 INTEGER INFO, KD, LDAB, LDAFB, LDB, LDX, N, NRHS
17 DOUBLE PRECISION RCOND
19 INTEGER IWORK( * )
21 DOUBLE PRECISION AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, *
23 ), BERR( * ), FERR( * ), S( * ), WORK( * ), X( LDX,
24 * )
25
27 DPBSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to com‐
28 pute the solution to a real system of linear equations
29 A * X = B, where A is an N-by-N symmetric positive definite band
30 matrix and X and B are N-by-NRHS matrices.
31 Error bounds on the solution and a condition estimate are also pro‐
33 vided.
34
37 The following steps are performed:
38 1. If FACT = 'E', real scaling factors are computed to equilibrate
40 the system:
41 diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * 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(S)*A*diag(S) and B by diag(S)*B.
45 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
47 factor the matrix A (after equilibration if FACT = 'E') as
48 A = U**T * U, if UPLO = 'U', or
49 A = L * L**T, if UPLO = 'L',
50 where U is an upper triangular band matrix, and L is a lower
51 triangular band matrix.
52 3. If the leading i-by-i principal minor is not positive definite,
54 then the routine returns with INFO = i. Otherwise, the factored
55 form of A is used to estimate the condition number of the matrix
56 A. If the reciprocal of the condition number is less than machine
57 precision, INFO = N+1 is returned as a warning, but the routine
58 still goes on to solve for X and compute error bounds as
59 described below.
60 4. The system of equations is solved for X using the factored form
62 of A.
63 5. Iterative refinement is applied to improve the computed solution
65 matrix and calculate error bounds and backward error estimates
66 for it.
67 6. If equilibration was used, the matrix X is premultiplied by
69 diag(S) so that it solves the original system before
70 equilibration.
71
74 FACT (input) CHARACTER*1
75 Specifies whether or not the factored form of the matrix A is
76 supplied on entry, and if not, whether the matrix A should be
77 equilibrated before it is factored. = 'F': On entry, AFB con‐
78 tains the factored form of A. If EQUED = 'Y', the matrix A has
79 been equilibrated with scaling factors given by S. AB and AFB
80 will not be modified. = 'N': The matrix A will be copied to
81 AFB and factored.
82 = 'E': The matrix A will be equilibrated if necessary, then
83 copied to AFB and factored.
84 UPLO (input) CHARACTER*1
86 = 'U': Upper triangle of A is stored;
87 = 'L': Lower triangle of A is stored.
88 N (input) INTEGER
90 The number of linear equations, i.e., the order of the matrix
91 A. N >= 0.
92 KD (input) INTEGER
94 The number of superdiagonals of the matrix A if UPLO = 'U', or
95 the number of subdiagonals if UPLO = 'L'. KD >= 0.
96 NRHS (input) INTEGER
98 The number of right-hand sides, i.e., the number of columns of
99 the matrices B and X. NRHS >= 0.
100 AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
102 On entry, the upper or lower triangle of the symmetric band
103 matrix A, stored in the first KD+1 rows of the array, except if
104 FACT = 'F' and EQUED = 'Y', then A must contain the equili‐
105 brated matrix diag(S)*A*diag(S). The j-th column of A is
106 stored in the j-th column of the array AB as follows: if UPLO =
107 'U', AB(KD+1+i-j,j) = A(i,j) for max(1,j-KD)<=i<=j; if UPLO =
108 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(N,j+KD). See below
109 for further details.
110 On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
112 diag(S)*A*diag(S).
113 LDAB (input) INTEGER
115 The leading dimension of the array A. LDAB >= KD+1.
116 AFB (input or output) DOUBLE PRECISION array, dimension (LDAFB,N)
118 If FACT = 'F', then AFB is an input argument and on entry con‐
119 tains the triangular factor U or L from the Cholesky factoriza‐
120 tion A = U**T*U or A = L*L**T of the band matrix A, in the same
121 storage format as A (see AB). If EQUED = 'Y', then AFB is the
122 factored form of the equilibrated matrix A.
123 If FACT = 'N', then AFB is an output argument and on exit
125 returns the triangular factor U or L from the Cholesky factor‐
126 ization A = U**T*U or A = L*L**T.
127 If FACT = 'E', then AFB is an output argument and on exit
129 returns the triangular factor U or L from the Cholesky factor‐
130 ization A = U**T*U or A = L*L**T of the equilibrated matrix A
131 (see the description of A for the form of the equilibrated
132 matrix).
133 LDAFB (input) INTEGER
135 The leading dimension of the array AFB. LDAFB >= KD+1.
136 EQUED (input or output) CHARACTER*1
138 Specifies the form of equilibration that was done. = 'N': No
139 equilibration (always true if FACT = 'N').
140 = 'Y': Equilibration was done, i.e., A has been replaced by
141 diag(S) * A * diag(S). EQUED is an input argument if FACT =
142 'F'; otherwise, it is an output argument.
143 S (input or output) DOUBLE PRECISION array, dimension (N)
145 The scale factors for A; not accessed if EQUED = 'N'. S is an
146 input argument if FACT = 'F'; otherwise, S is an output argu‐
147 ment. If FACT = 'F' and EQUED = 'Y', each element of S must be
148 positive.
149 B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
151 On entry, the N-by-NRHS right hand side matrix B. On exit, if
152 EQUED = 'N', B is not modified; if EQUED = 'Y', B is overwrit‐
153 ten by diag(S) * B.
154 LDB (input) INTEGER
156 The leading dimension of the array B. LDB >= max(1,N).
157 X (output) DOUBLE PRECISION array, dimension (LDX,NRHS)
159 If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
160 the original system of equations. Note that if EQUED = 'Y', A
161 and B are modified on exit, and the solution to the equili‐
162 brated system is inv(diag(S))*X.
163 LDX (input) INTEGER
165 The leading dimension of the array X. LDX >= max(1,N).
166 RCOND (output) DOUBLE PRECISION
168 The estimate of the reciprocal condition number of the matrix A
169 after equilibration (if done). If RCOND is less than the
170 machine precision (in particular, if RCOND = 0), the matrix is
171 singular to working precision. This condition is indicated by
172 a return code of INFO > 0.
173 FERR (output) DOUBLE PRECISION array, dimension (NRHS)
175 The estimated forward error bound for each solution vector X(j)
176 (the j-th column of the solution matrix X). If XTRUE is the
177 true solution corresponding to X(j), FERR(j) is an estimated
178 upper bound for the magnitude of the largest element in (X(j) -
179 XTRUE) divided by the magnitude of the largest element in X(j).
180 The estimate is as reliable as the estimate for RCOND, and is
181 almost always a slight overestimate of the true error.
182 BERR (output) DOUBLE PRECISION array, dimension (NRHS)
184 The componentwise relative backward error of each solution vec‐
185 tor X(j) (i.e., the smallest relative change in any element of
186 A or B that makes X(j) an exact solution).
187 WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
189 IWORK (workspace) INTEGER array, dimension (N)
191 INFO (output) INTEGER
193 = 0: successful exit
194 < 0: if INFO = -i, the i-th argument had an illegal value
195 > 0: if INFO = i, and i is
196 <= N: the leading minor of order i of A is not positive defi‐
197 nite, so the factorization could not be completed, and the
198 solution has not been computed. RCOND = 0 is returned. = N+1:
199 U is nonsingular, but RCOND is less than machine precision,
200 meaning that the matrix is singular to working precision. Nev‐
201 ertheless, the solution and error bounds are computed because
202 there are a number of situations where the computed solution
203 can be more accurate than the value of RCOND would suggest.
204
206 The band storage scheme is illustrated by the following example, when N
207 = 6, KD = 2, and UPLO = 'U':
208 Two-dimensional storage of the symmetric matrix A:
210 a11 a12 a13
212 a22 a23 a24
213 a33 a34 a35
214 a44 a45 a46
215 a55 a56
216 (aij=conjg(aji)) a66
217 Band storage of the upper triangle of A:
219 * * a13 a24 a35 a46
221 * a12 a23 a34 a45 a56
222 a11 a22 a33 a44 a55 a66
223 Similarly, if UPLO = 'L' the format of A is as follows:
225 a11 a22 a33 a44 a55 a66
227 a21 a32 a43 a54 a65 *
228 a31 a42 a53 a64 * *
229 Array elements marked * are not used by the routine.
231 LAPACK driver routine (version 3.N1o)vember 2006 DPBSVX(1)