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