1DPOSVX(1) LAPACK driver routine (version 3.2) DPOSVX(1)
2
6 DPOSVX - 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 DPOSVX( 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 DOUBLE PRECISION RCOND
18 INTEGER IWORK( * )
20 DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
22 BERR( * ), FERR( * ), S( * ), WORK( * ), X( LDX, * )
23
25 DPOSVX 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‐
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, AF con‐
67 tains the factored form of A. If EQUED = 'Y', the matrix A has
68 been equilibrated with scaling factors given by S. A and AF
69 will not be modified. = 'N': The matrix A will be copied to
70 AF and factored.
71 = 'E': The matrix A will be equilibrated if necessary, then
72 copied to AF 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 A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
87 On entry, the symmetric matrix A, except if FACT = 'F' and
88 EQUED = 'Y', then A must contain the equilibrated matrix
89 diag(S)*A*diag(S). If UPLO = 'U', the leading N-by-N upper
90 triangular part of A contains the upper triangular part of the
91 matrix A, and the strictly lower triangular part of A is not
92 referenced. If UPLO = 'L', the leading N-by-N lower triangular
93 part of A contains the lower triangular part of the matrix A,
94 and the strictly upper triangular part of A is not referenced.
95 A is not modified if FACT = 'F' or 'N', or if FACT = 'E' and
96 EQUED = 'N' on exit. On exit, if FACT = 'E' and EQUED = 'Y', A
97 is overwritten by diag(S)*A*diag(S).
98 LDA (input) INTEGER
100 The leading dimension of the array A. LDA >= max(1,N).
101 AF (input or output) DOUBLE PRECISION array, dimension (LDAF,N)
103 If FACT = 'F', then AF is an input argument and on entry con‐
104 tains the triangular factor U or L from the Cholesky factoriza‐
105 tion A = U**T*U or A = L*L**T, in the same storage format as A.
106 If EQUED .ne. 'N', then AF is the factored form of the equili‐
107 brated matrix diag(S)*A*diag(S). If FACT = 'N', then AF is an
108 output argument and on exit returns the triangular factor U or
109 L from the Cholesky factorization A = U**T*U or A = L*L**T of
110 the original matrix A. If FACT = 'E', then AF is an output
111 argument and on exit returns the triangular factor U or L from
112 the Cholesky factorization A = U**T*U or A = L*L**T of the
113 equilibrated matrix A (see the description of A for the form of
114 the equilibrated matrix).
115 LDAF (input) INTEGER
117 The leading dimension of the array AF. LDAF >= max(1,N).
118 EQUED (input or output) CHARACTER*1
120 Specifies the form of equilibration that was done. = 'N': No
121 equilibration (always true if FACT = 'N').
122 = 'Y': Equilibration was done, i.e., A has been replaced by
123 diag(S) * A * diag(S). EQUED is an input argument if FACT =
124 'F'; otherwise, it is an output argument.
125 S (input or output) DOUBLE PRECISION array, dimension (N)
127 The scale factors for A; not accessed if EQUED = 'N'. S is an
128 input argument if FACT = 'F'; otherwise, S is an output argu‐
129 ment. If FACT = 'F' and EQUED = 'Y', each element of S must be
130 positive.
131 B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
133 On entry, the N-by-NRHS right hand side matrix B. On exit, if
134 EQUED = 'N', B is not modified; if EQUED = 'Y', B is overwrit‐
135 ten by diag(S) * B.
136 LDB (input) INTEGER
138 The leading dimension of the array B. LDB >= max(1,N).
139 X (output) DOUBLE PRECISION array, dimension (LDX,NRHS)
141 If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
142 the original system of equations. Note that if EQUED = 'Y', A
143 and B are modified on exit, and the solution to the equili‐
144 brated system is inv(diag(S))*X.
145 LDX (input) INTEGER
147 The leading dimension of the array X. LDX >= max(1,N).
148 RCOND (output) DOUBLE PRECISION
150 The estimate of the reciprocal condition number of the matrix A
151 after equilibration (if done). If RCOND is less than the
152 machine precision (in particular, if RCOND = 0), the matrix is
153 singular to working precision. This condition is indicated by
154 a return code of INFO > 0.
155 FERR (output) DOUBLE PRECISION array, dimension (NRHS)
157 The estimated forward error bound for each solution vector X(j)
158 (the j-th column of the solution matrix X). If XTRUE is the
159 true solution corresponding to X(j), FERR(j) is an estimated
160 upper bound for the magnitude of the largest element in (X(j) -
161 XTRUE) divided by the magnitude of the largest element in X(j).
162 The estimate is as reliable as the estimate for RCOND, and is
163 almost always a slight overestimate of the true error.
164 BERR (output) DOUBLE PRECISION array, dimension (NRHS)
166 The componentwise relative backward error of each solution vec‐
167 tor X(j) (i.e., the smallest relative change in any element of
168 A or B that makes X(j) an exact solution).
169 WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
171 IWORK (workspace) INTEGER array, dimension (N)
173 INFO (output) INTEGER
175 = 0: successful exit
176 < 0: if INFO = -i, the i-th argument had an illegal value
177 > 0: if INFO = i, and i is
178 <= N: the leading minor of order i of A is not positive defi‐
179 nite, so the factorization could not be completed, and the
180 solution has not been computed. RCOND = 0 is returned. = N+1:
181 U is nonsingular, but RCOND is less than machine precision,
182 meaning that the matrix is singular to working precision. Nev‐
183 ertheless, the solution and error bounds are computed because
184 there are a number of situations where the computed solution
185 can be more accurate than the value of RCOND would suggest.
186 LAPACK driver routine (version 3.N2o)vember 2008 DPOSVX(1)