1SPOSVX(1) LAPACK driver routine (version 3.1) SPOSVX(1)
2
6 SPOSVX - 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 SPOSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED, S, B,
11 LDB, X, LDX, RCOND, FERR, BERR, WORK, IWORK, INFO )
12 CHARACTER EQUED, FACT, UPLO
14 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS
16 REAL RCOND
18 INTEGER IWORK( * )
20 REAL A( LDA, * ), AF( LDAF, * ), B( LDB, * ), BERR( * ),
22 FERR( * ), S( * ), WORK( * ), X( LDX, * )
23
25 SPOSVX 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 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, AF con‐
76 tains the factored form of A. If EQUED = 'Y', the matrix A has
77 been equilibrated with scaling factors given by S. A and AF
78 will not be modified. = 'N': The matrix A will be copied to
79 AF and factored.
80 = 'E': The matrix A will be equilibrated if necessary, then
81 copied to AF 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 A (input/output) REAL array, dimension (LDA,N)
96 On entry, the symmetric matrix A, except if FACT = 'F' and
97 EQUED = 'Y', then A must contain the equilibrated matrix
98 diag(S)*A*diag(S). If UPLO = 'U', the leading N-by-N upper
99 triangular part of A contains the upper triangular part of the
100 matrix A, and the strictly lower triangular part of A is not
101 referenced. If UPLO = 'L', the leading N-by-N lower triangular
102 part of A contains the lower triangular part of the matrix A,
103 and the strictly upper triangular part of A is not referenced.
104 A is not modified if FACT = 'F' or 'N', or if FACT = 'E' and
105 EQUED = 'N' on exit.
106 On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
108 diag(S)*A*diag(S).
109 LDA (input) INTEGER
111 The leading dimension of the array A. LDA >= max(1,N).
112 AF (input or output) REAL array, dimension (LDAF,N)
114 If FACT = 'F', then AF is an input argument and on entry con‐
115 tains the triangular factor U or L from the Cholesky factoriza‐
116 tion A = U**T*U or A = L*L**T, in the same storage format as A.
117 If EQUED .ne. 'N', then AF is the factored form of the equili‐
118 brated matrix diag(S)*A*diag(S).
119 If FACT = 'N', then AF is an output argument and on exit
121 returns the triangular factor U or L from the Cholesky factor‐
122 ization A = U**T*U or A = L*L**T of the original matrix A.
123 If FACT = 'E', then AF is an output argument and on exit
125 returns the triangular factor U or L from the Cholesky factor‐
126 ization A = U**T*U or A = L*L**T of the equilibrated matrix A
127 (see the description of A for the form of the equilibrated
128 matrix).
129 LDAF (input) INTEGER
131 The leading dimension of the array AF. LDAF >= max(1,N).
132 EQUED (input or output) CHARACTER*1
134 Specifies the form of equilibration that was done. = 'N': No
135 equilibration (always true if FACT = 'N').
136 = 'Y': Equilibration was done, i.e., A has been replaced by
137 diag(S) * A * diag(S). EQUED is an input argument if FACT =
138 'F'; otherwise, it is an output argument.
139 S (input or output) REAL array, dimension (N)
141 The scale factors for A; not accessed if EQUED = 'N'. S is an
142 input argument if FACT = 'F'; otherwise, S is an output argu‐
143 ment. If FACT = 'F' and EQUED = 'Y', each element of S must be
144 positive.
145 B (input/output) REAL array, dimension (LDB,NRHS)
147 On entry, the N-by-NRHS right hand side matrix B. On exit, if
148 EQUED = 'N', B is not modified; if EQUED = 'Y', B is overwrit‐
149 ten by diag(S) * B.
150 LDB (input) INTEGER
152 The leading dimension of the array B. LDB >= max(1,N).
153 X (output) REAL array, dimension (LDX,NRHS)
155 If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
156 the original system of equations. Note that if EQUED = 'Y', A
157 and B are modified on exit, and the solution to the equili‐
158 brated system is inv(diag(S))*X.
159 LDX (input) INTEGER
161 The leading dimension of the array X. LDX >= max(1,N).
162 RCOND (output) REAL
164 The estimate of the reciprocal condition number of the matrix A
165 after equilibration (if done). If RCOND is less than the
166 machine precision (in particular, if RCOND = 0), the matrix is
167 singular to working precision. This condition is indicated by
168 a return code of INFO > 0.
169 FERR (output) REAL array, dimension (NRHS)
171 The estimated forward error bound for each solution vector X(j)
172 (the j-th column of the solution matrix X). If XTRUE is the
173 true solution corresponding to X(j), FERR(j) is an estimated
174 upper bound for the magnitude of the largest element in (X(j) -
175 XTRUE) divided by the magnitude of the largest element in X(j).
176 The estimate is as reliable as the estimate for RCOND, and is
177 almost always a slight overestimate of the true error.
178 BERR (output) REAL array, dimension (NRHS)
180 The componentwise relative backward error of each solution vec‐
181 tor X(j) (i.e., the smallest relative change in any element of
182 A or B that makes X(j) an exact solution).
183 WORK (workspace) REAL array, dimension (3*N)
185 IWORK (workspace) INTEGER array, dimension (N)
187 INFO (output) INTEGER
189 = 0: successful exit
190 < 0: if INFO = -i, the i-th argument had an illegal value
191 > 0: if INFO = i, and i is
192 <= N: the leading minor of order i of A is not positive defi‐
193 nite, so the factorization could not be completed, and the
194 solution has not been computed. RCOND = 0 is returned. = N+1:
195 U is nonsingular, but RCOND is less than machine precision,
196 meaning that the matrix is singular to working precision. Nev‐
197 ertheless, the solution and error bounds are computed because
198 there are a number of situations where the computed solution
199 can be more accurate than the value of RCOND would suggest.
200 LAPACK driver routine (version 3.N1o)vember 2006 SPOSVX(1)