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