1CPPSVX(1) LAPACK driver routine (version 3.1) CPPSVX(1)
2
6 CPPSVX - the Cholesky factorization A = U**H*U or A = L*L**H to compute
7 the solution to a complex system of linear equations A * X = B,
8
10 SUBROUTINE CPPSVX( FACT, UPLO, N, NRHS, AP, AFP, EQUED, S, B, LDB, X,
11 LDX, RCOND, FERR, BERR, WORK, RWORK, INFO )
12 CHARACTER EQUED, FACT, UPLO
14 INTEGER INFO, LDB, LDX, N, NRHS
16 REAL RCOND
18 REAL BERR( * ), FERR( * ), RWORK( * ), S( * )
20 COMPLEX AFP( * ), AP( * ), B( LDB, * ), WORK( * ), X( LDX, *
22 )
23
25 CPPSVX uses the Cholesky factorization A = U**H*U or A = L*L**H to com‐
26 pute the solution to a complex system of linear equations
27 A * X = B, where A is an N-by-N Hermitian 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'* U , if UPLO = 'U', or
47 A = L * L', if UPLO = 'L',
48 where U is an upper triangular matrix, L is a lower triangular
49 matrix, and ' indicates conjugate transpose.
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) COMPLEX array, dimension (N*(N+1)/2)
96 On entry, the upper or lower triangle of the Hermitian 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) COMPLEX array, dimension (N*(N+1)/2)
109 If FACT = 'F', then AFP is an input argument and on entry con‐
110 tains the triangular factor U or L from the Cholesky factoriza‐
111 tion A = U**H*U or A = L*L**H, in the same storage format as A.
112 If EQUED .ne. 'N', then AFP is the factored form of the equili‐
113 brated 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**H*U or A = L*L**H 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**H*U or A = L*L**H of the equilibrated matrix A
122 (see the description of AP for the form of the equilibrated
123 matrix).
124 EQUED (input or output) CHARACTER*1
126 Specifies the form of equilibration that was done. = 'N': No
127 equilibration (always true if FACT = 'N').
128 = 'Y': Equilibration was done, i.e., A has been replaced by
129 diag(S) * A * diag(S). EQUED is an input argument if FACT =
130 'F'; otherwise, it is an output argument.
131 S (input or output) REAL array, dimension (N)
133 The scale factors for A; not accessed if EQUED = 'N'. S is an
134 input argument if FACT = 'F'; otherwise, S is an output argu‐
135 ment. If FACT = 'F' and EQUED = 'Y', each element of S must be
136 positive.
137 B (input/output) COMPLEX array, dimension (LDB,NRHS)
139 On entry, the N-by-NRHS right hand side matrix B. On exit, if
140 EQUED = 'N', B is not modified; if EQUED = 'Y', B is overwrit‐
141 ten by diag(S) * B.
142 LDB (input) INTEGER
144 The leading dimension of the array B. LDB >= max(1,N).
145 X (output) COMPLEX array, dimension (LDX,NRHS)
147 If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
148 the original system of equations. Note that if EQUED = 'Y', A
149 and B are modified on exit, and the solution to the equili‐
150 brated system is inv(diag(S))*X.
151 LDX (input) INTEGER
153 The leading dimension of the array X. LDX >= max(1,N).
154 RCOND (output) REAL
156 The estimate of the reciprocal condition number of the matrix A
157 after equilibration (if done). If RCOND is less than the
158 machine precision (in particular, if RCOND = 0), the matrix is
159 singular to working precision. This condition is indicated by
160 a return code of INFO > 0.
161 FERR (output) REAL array, dimension (NRHS)
163 The estimated forward error bound for each solution vector X(j)
164 (the j-th column of the solution matrix X). If XTRUE is the
165 true solution corresponding to X(j), FERR(j) is an estimated
166 upper bound for the magnitude of the largest element in (X(j) -
167 XTRUE) divided by the magnitude of the largest element in X(j).
168 The estimate is as reliable as the estimate for RCOND, and is
169 almost always a slight overestimate of the true error.
170 BERR (output) REAL array, dimension (NRHS)
172 The componentwise relative backward error of each solution vec‐
173 tor X(j) (i.e., the smallest relative change in any element of
174 A or B that makes X(j) an exact solution).
175 WORK (workspace) COMPLEX array, dimension (2*N)
177 RWORK (workspace) REAL array, dimension (N)
179 INFO (output) INTEGER
181 = 0: successful exit
182 < 0: if INFO = -i, the i-th argument had an illegal value
183 > 0: if INFO = i, and i is
184 <= N: the leading minor of order i of A is not positive defi‐
185 nite, so the factorization could not be completed, and the
186 solution has not been computed. RCOND = 0 is returned. = N+1:
187 U is nonsingular, but RCOND is less than machine precision,
188 meaning that the matrix is singular to working precision. Nev‐
189 ertheless, the solution and error bounds are computed because
190 there are a number of situations where the computed solution
191 can be more accurate than the value of RCOND would suggest.
192
194 The packed storage scheme is illustrated by the following example when
195 N = 4, UPLO = 'U':
196 Two-dimensional storage of the Hermitian matrix A:
198 a11 a12 a13 a14
200 a22 a23 a24
201 a33 a34 (aij = conjg(aji))
202 a44
203 Packed storage of the upper triangle of A:
205 AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
207 LAPACK driver routine (version 3.N1o)vember 2006 CPPSVX(1)