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