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