1SPPSVX(1) LAPACK driver routine (version 3.1) SPPSVX(1)
2
6 SPPSVX - 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 SPPSVX( FACT, UPLO, N, NRHS, AP, AFP, EQUED, S, B, LDB, X,
11 LDX, RCOND, FERR, BERR, WORK, IWORK, INFO )
12 CHARACTER EQUED, FACT, UPLO
14 INTEGER INFO, LDB, LDX, N, NRHS
16 REAL RCOND
18 INTEGER IWORK( * )
20 REAL AFP( * ), AP( * ), B( LDB, * ), BERR( * ), FERR( *
22 ), S( * ), WORK( * ), X( LDX, * )
23
25 SPPSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to com‐
26 pute the solution to a real system of linear equations
27 A * X = B, where A is an N-by-N symmetric positive definite matrix
28 stored in packed format and X and B are N-by-NRHS matrices.
29 Error bounds on the solution and a condition estimate are also pro‐
31 vided.
32
35 The following steps are performed:
36 1. If FACT = 'E', real scaling factors are computed to equilibrate
38 the system:
39 diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B
40 Whether or not the system will be equilibrated depends on the
41 scaling of the matrix A, but if equilibration is used, A is
42 overwritten by diag(S)*A*diag(S) and B by diag(S)*B.
43 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
45 factor the matrix A (after equilibration if FACT = 'E') as
46 A = U**T* U, if UPLO = 'U', or
47 A = L * L**T, if UPLO = 'L',
48 where U is an upper triangular matrix and L is a lower triangular
49 matrix.
50 3. If the leading i-by-i principal minor is not positive definite,
52 then the routine returns with INFO = i. Otherwise, the factored
53 form of A is used to estimate the condition number of the matrix
54 A. If the reciprocal of the condition number is less than machine
55 precision, INFO = N+1 is returned as a warning, but the routine
56 still goes on to solve for X and compute error bounds as
57 described below.
58 4. The 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
63 matrix and calculate error bounds and backward error estimates
64 for it.
65 6. If equilibration was used, the matrix X is premultiplied by
67 diag(S) so that it solves the original system before
68 equilibration.
69
72 FACT (input) CHARACTER*1
73 Specifies whether or not the factored form of the matrix A is
74 supplied on entry, and if not, whether the matrix A should be
75 equilibrated before it is factored. = 'F': On entry, AFP con‐
76 tains the factored form of A. If EQUED = 'Y', the matrix A has
77 been equilibrated with scaling factors given by S. AP and AFP
78 will not be modified. = 'N': The matrix A will be copied to
79 AFP and factored.
80 = 'E': The matrix A will be equilibrated if necessary, then
81 copied to AFP and factored.
82 UPLO (input) CHARACTER*1
84 = 'U': Upper triangle of A is stored;
85 = 'L': Lower triangle of A is stored.
86 N (input) INTEGER
88 The number of linear equations, i.e., the order of the matrix
89 A. N >= 0.
90 NRHS (input) INTEGER
92 The number of right hand sides, i.e., the number of columns of
93 the matrices B and X. NRHS >= 0.
94 AP (input/output) REAL array, dimension (N*(N+1)/2)
96 On entry, the upper or lower triangle of the symmetric matrix
97 A, packed columnwise in a linear array, except if FACT = 'F'
98 and EQUED = 'Y', then A must contain the equilibrated matrix
99 diag(S)*A*diag(S). The j-th column of A is stored in the array
100 AP as follows: if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for
101 1<=i<=j; if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for
102 j<=i<=n. See below for further details. A is not modified if
103 FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit.
104 On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
106 diag(S)*A*diag(S).
107 AFP (input or output) REAL array, dimension
109 (N*(N+1)/2) If FACT = 'F', then AFP is an input argument and on
110 entry contains the triangular factor U or L from the Cholesky
111 factorization A = U'*U or A = L*L', in the same storage format
112 as A. If EQUED .ne. 'N', then AFP is the factored form of the
113 equilibrated matrix A.
114 If FACT = 'N', then AFP is an output argument and on exit
116 returns the triangular factor U or L from the Cholesky factor‐
117 ization A = U'*U or A = L*L' of the original matrix A.
118 If FACT = 'E', then AFP is an output argument and on exit
120 returns the triangular factor U or L from the Cholesky factor‐
121 ization A = U'*U or A = L*L' of the equilibrated matrix A (see
122 the description of AP for the form of the equilibrated matrix).
123 EQUED (input or output) CHARACTER*1
125 Specifies the form of equilibration that was done. = 'N': No
126 equilibration (always true if FACT = 'N').
127 = 'Y': Equilibration was done, i.e., A has been replaced by
128 diag(S) * A * diag(S). EQUED is an input argument if FACT =
129 'F'; otherwise, it is an output argument.
130 S (input or output) REAL array, dimension (N)
132 The scale factors for A; not accessed if EQUED = 'N'. S is an
133 input argument if FACT = 'F'; otherwise, S is an output argu‐
134 ment. If FACT = 'F' and EQUED = 'Y', each element of S must be
135 positive.
136 B (input/output) REAL array, dimension (LDB,NRHS)
138 On entry, the N-by-NRHS right hand side matrix B. On exit, if
139 EQUED = 'N', B is not modified; if EQUED = 'Y', B is overwrit‐
140 ten by diag(S) * B.
141 LDB (input) INTEGER
143 The leading dimension of the array B. LDB >= max(1,N).
144 X (output) REAL array, dimension (LDX,NRHS)
146 If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
147 the original system of equations. Note that if EQUED = 'Y', A
148 and B are modified on exit, and the solution to the equili‐
149 brated system is inv(diag(S))*X.
150 LDX (input) INTEGER
152 The leading dimension of the array X. LDX >= max(1,N).
153 RCOND (output) REAL
155 The estimate of the reciprocal condition number of the matrix A
156 after equilibration (if done). If RCOND is less than the
157 machine precision (in particular, if RCOND = 0), the matrix is
158 singular to working precision. This condition is indicated by
159 a return code of INFO > 0.
160 FERR (output) REAL array, dimension (NRHS)
162 The estimated forward error bound for each solution vector X(j)
163 (the j-th column of the solution matrix X). If XTRUE is the
164 true solution corresponding to X(j), FERR(j) is an estimated
165 upper bound for the magnitude of the largest element in (X(j) -
166 XTRUE) divided by the magnitude of the largest element in X(j).
167 The estimate is as reliable as the estimate for RCOND, and is
168 almost always a slight overestimate of the true error.
169 BERR (output) REAL array, dimension (NRHS)
171 The componentwise relative backward error of each solution vec‐
172 tor X(j) (i.e., the smallest relative change in any element of
173 A or B that makes X(j) an exact solution).
174 WORK (workspace) REAL array, dimension (3*N)
176 IWORK (workspace) INTEGER array, dimension (N)
178 INFO (output) INTEGER
180 = 0: successful exit
181 < 0: if INFO = -i, the i-th argument had an illegal value
182 > 0: if INFO = i, and i is
183 <= N: the leading minor of order i of A is not positive defi‐
184 nite, so the factorization could not be completed, and the
185 solution has not been computed. RCOND = 0 is returned. = N+1:
186 U is nonsingular, but RCOND is less than machine precision,
187 meaning that the matrix is singular to working precision. Nev‐
188 ertheless, the solution and error bounds are computed because
189 there are a number of situations where the computed solution
190 can be more accurate than the value of RCOND would suggest.
191
193 The packed storage scheme is illustrated by the following example when
194 N = 4, UPLO = 'U':
195 Two-dimensional storage of the symmetric matrix A:
197 a11 a12 a13 a14
199 a22 a23 a24
200 a33 a34 (aij = conjg(aji))
201 a44
202 Packed storage of the upper triangle of A:
204 AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
206 LAPACK driver routine (version 3.N1o)vember 2006 SPPSVX(1)