1ZPOSVX(1) LAPACK driver routine (version 3.2) ZPOSVX(1)
2
6 ZPOSVX - uses the Cholesky factorization A = U**H*U or A = L*L**H to
7 compute the solution to a complex system of linear equations A * X =
8 B,
9
11 SUBROUTINE ZPOSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED, S, B,
12 LDB, X, LDX, RCOND, FERR, BERR, WORK, RWORK, INFO )
13 CHARACTER EQUED, FACT, UPLO
15 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS
17 DOUBLE PRECISION RCOND
19 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * ), S( * )
21 COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ), WORK( * ),
23 X( LDX, * )
24
26 ZPOSVX uses the Cholesky factorization A = U**H*U or A = L*L**H to com‐
27 pute the solution to a complex system of linear equations
28 A * X = B, where A is an N-by-N Hermitian positive definite matrix
29 and X and B are N-by-NRHS matrices.
30 Error bounds on the solution and a condition estimate are also pro‐
31 vided.
32
34 The following steps are performed:
35 1. If FACT = 'E', real scaling factors are computed to equilibrate
36 the system:
37 diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B
38 Whether or not the system will be equilibrated depends on the
39 scaling of the matrix A, but if equilibration is used, A is
40 overwritten by diag(S)*A*diag(S) and B by diag(S)*B.
41 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
42 factor the matrix A (after equilibration if FACT = 'E') as
43 A = U**H* U, if UPLO = 'U', or
44 A = L * L**H, if UPLO = 'L',
45 where U is an upper triangular matrix and L is a lower triangular
46 matrix.
47 3. If the leading i-by-i principal minor is not positive definite,
48 then the routine returns with INFO = i. Otherwise, the factored
49 form of A is used to estimate the condition number of the matrix
50 A. If the reciprocal of the condition number is less than machine
51 precision, INFO = N+1 is returned as a warning, but the routine
52 still goes on to solve for X and compute error bounds as
53 described below.
54 4. The system of equations is solved for X using the factored form
55 of A.
56 5. Iterative refinement is applied to improve the computed solution
57 matrix and calculate error bounds and backward error estimates
58 for it.
59 6. If equilibration was used, the matrix X is premultiplied by
60 diag(S) so that it solves the original system before
61 equilibration.
62
64 FACT (input) CHARACTER*1
65 Specifies whether or not the factored form of the matrix A is
66 supplied on entry, and if not, whether the matrix A should be
67 equilibrated before it is factored. = 'F': On entry, AF con‐
68 tains the factored form of A. If EQUED = 'Y', the matrix A has
69 been equilibrated with scaling factors given by S. A and AF
70 will not be modified. = 'N': The matrix A will be copied to
71 AF and factored.
72 = 'E': The matrix A will be equilibrated if necessary, then
73 copied to AF and factored.
74 UPLO (input) CHARACTER*1
76 = 'U': Upper triangle of A is stored;
77 = 'L': Lower triangle of A is stored.
78 N (input) INTEGER
80 The number of linear equations, i.e., the order of the matrix
81 A. N >= 0.
82 NRHS (input) INTEGER
84 The number of right hand sides, i.e., the number of columns of
85 the matrices B and X. NRHS >= 0.
86 A (input/output) COMPLEX*16 array, dimension (LDA,N)
88 On entry, the Hermitian matrix A, except if FACT = 'F' and
89 EQUED = 'Y', then A must contain the equilibrated matrix
90 diag(S)*A*diag(S). If UPLO = 'U', the leading N-by-N upper
91 triangular part of A contains the upper triangular part of the
92 matrix A, and the strictly lower triangular part of A is not
93 referenced. If UPLO = 'L', the leading N-by-N lower triangular
94 part of A contains the lower triangular part of the matrix A,
95 and the strictly upper triangular part of A is not referenced.
96 A is not modified if FACT = 'F' or 'N', or if FACT = 'E' and
97 EQUED = 'N' on exit. On exit, if FACT = 'E' and EQUED = 'Y', A
98 is overwritten by diag(S)*A*diag(S).
99 LDA (input) INTEGER
101 The leading dimension of the array A. LDA >= max(1,N).
102 AF (input or output) COMPLEX*16 array, dimension (LDAF,N)
104 If FACT = 'F', then AF is an input argument and on entry con‐
105 tains the triangular factor U or L from the Cholesky factoriza‐
106 tion A = U**H*U or A = L*L**H, in the same storage format as A.
107 If EQUED .ne. 'N', then AF is the factored form of the equili‐
108 brated matrix diag(S)*A*diag(S). If FACT = 'N', then AF is an
109 output argument and on exit returns the triangular factor U or
110 L from the Cholesky factorization A = U**H*U or A = L*L**H of
111 the original matrix A. If FACT = 'E', then AF is an output
112 argument and on exit returns the triangular factor U or L from
113 the Cholesky factorization A = U**H*U or A = L*L**H of the
114 equilibrated matrix A (see the description of A for the form of
115 the equilibrated matrix).
116 LDAF (input) INTEGER
118 The leading dimension of the array AF. LDAF >= max(1,N).
119 EQUED (input or output) CHARACTER*1
121 Specifies the form of equilibration that was done. = 'N': No
122 equilibration (always true if FACT = 'N').
123 = 'Y': Equilibration was done, i.e., A has been replaced by
124 diag(S) * A * diag(S). EQUED is an input argument if FACT =
125 'F'; otherwise, it is an output argument.
126 S (input or output) DOUBLE PRECISION array, dimension (N)
128 The scale factors for A; not accessed if EQUED = 'N'. S is an
129 input argument if FACT = 'F'; otherwise, S is an output argu‐
130 ment. If FACT = 'F' and EQUED = 'Y', each element of S must be
131 positive.
132 B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
134 On entry, the N-by-NRHS righthand side matrix B. On exit, if
135 EQUED = 'N', B is not modified; if EQUED = 'Y', B is overwrit‐
136 ten by diag(S) * B.
137 LDB (input) INTEGER
139 The leading dimension of the array B. LDB >= max(1,N).
140 X (output) COMPLEX*16 array, dimension (LDX,NRHS)
142 If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
143 the original system of equations. Note that if EQUED = 'Y', A
144 and B are modified on exit, and the solution to the equili‐
145 brated system is inv(diag(S))*X.
146 LDX (input) INTEGER
148 The leading dimension of the array X. LDX >= max(1,N).
149 RCOND (output) DOUBLE PRECISION
151 The estimate of the reciprocal condition number of the matrix A
152 after equilibration (if done). If RCOND is less than the
153 machine precision (in particular, if RCOND = 0), the matrix is
154 singular to working precision. This condition is indicated by
155 a return code of INFO > 0.
156 FERR (output) DOUBLE PRECISION array, dimension (NRHS)
158 The estimated forward error bound for each solution vector X(j)
159 (the j-th column of the solution matrix X). If XTRUE is the
160 true solution corresponding to X(j), FERR(j) is an estimated
161 upper bound for the magnitude of the largest element in (X(j) -
162 XTRUE) divided by the magnitude of the largest element in X(j).
163 The estimate is as reliable as the estimate for RCOND, and is
164 almost always a slight overestimate of the true error.
165 BERR (output) DOUBLE PRECISION array, dimension (NRHS)
167 The componentwise relative backward error of each solution vec‐
168 tor X(j) (i.e., the smallest relative change in any element of
169 A or B that makes X(j) an exact solution).
170 WORK (workspace) COMPLEX*16 array, dimension (2*N)
172 RWORK (workspace) DOUBLE PRECISION array, dimension (N)
174 INFO (output) INTEGER
176 = 0: successful exit
177 < 0: if INFO = -i, the i-th argument had an illegal value
178 > 0: if INFO = i, and i is
179 <= N: the leading minor of order i of A is not positive defi‐
180 nite, so the factorization could not be completed, and the
181 solution has not been computed. RCOND = 0 is returned. = N+1:
182 U is nonsingular, but RCOND is less than machine precision,
183 meaning that the matrix is singular to working precision. Nev‐
184 ertheless, the solution and error bounds are computed because
185 there are a number of situations where the computed solution
186 can be more accurate than the value of RCOND would suggest.
187 LAPACK driver routine (version 3.N2o)vember 2008 ZPOSVX(1)