1ZSPSVX(1) LAPACK driver routine (version 3.1) ZSPSVX(1)
2
6 ZSPSVX - the diagonal pivoting factorization A = U*D*U**T or A =
7 L*D*L**T to compute the solution to a complex system of linear equa‐
8 tions A * X = B, where A is an N-by-N symmetric matrix stored in packed
9 format and X and B are N-by-NRHS matrices
10
12 SUBROUTINE ZSPSVX( FACT, UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X, LDX,
13 RCOND, FERR, BERR, WORK, RWORK, INFO )
14 CHARACTER FACT, UPLO
16 INTEGER INFO, LDB, LDX, N, NRHS
18 DOUBLE PRECISION RCOND
20 INTEGER IPIV( * )
22 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
24 COMPLEX*16 AFP( * ), AP( * ), B( LDB, * ), WORK( * ), X( LDX, *
26 )
27
29 ZSPSVX uses the diagonal pivoting factorization A = U*D*U**T or A =
30 L*D*L**T to compute the solution to a complex system of linear equa‐
31 tions A * X = B, where A is an N-by-N symmetric matrix stored in packed
32 format and X and B are N-by-NRHS matrices.
33 Error bounds on the solution and a condition estimate are also pro‐
35 vided.
36
39 The following steps are performed:
40 1. If FACT = 'N', the diagonal pivoting method is used to factor A as
42 A = U * D * U**T, if UPLO = 'U', or
43 A = L * D * L**T, if UPLO = 'L',
44 where U (or L) is a product of permutation and unit upper (lower)
45 triangular matrices and D is symmetric and block diagonal with
46 1-by-1 and 2-by-2 diagonal blocks.
47 2. If some D(i,i)=0, so that D is exactly singular, then the routine
49 returns with INFO = i. Otherwise, the factored form of A is used
50 to estimate the condition number of the matrix A. If the
51 reciprocal of the condition number is less than machine precision,
52 INFO = N+1 is returned as a warning, but the routine still goes on
53 to solve for X and compute error bounds as described below.
54 3. The system of equations is solved for X using the factored form
56 of A.
57 4. Iterative refinement is applied to improve the computed solution
59 matrix and calculate error bounds and backward error estimates
60 for it.
61
64 FACT (input) CHARACTER*1
65 Specifies whether or not the factored form of A has been sup‐
66 plied on entry. = 'F': On entry, AFP and IPIV contain the
67 factored form of A. AP, AFP and IPIV will not be modified. =
68 'N': The matrix A will be copied to AFP and factored.
69 UPLO (input) CHARACTER*1
71 = 'U': Upper triangle of A is stored;
72 = 'L': Lower triangle of A is stored.
73 N (input) INTEGER
75 The number of linear equations, i.e., the order of the matrix
76 A. N >= 0.
77 NRHS (input) INTEGER
79 The number of right hand sides, i.e., the number of columns of
80 the matrices B and X. NRHS >= 0.
81 AP (input) COMPLEX*16 array, dimension (N*(N+1)/2)
83 The upper or lower triangle of the symmetric matrix A, packed
84 columnwise in a linear array. The j-th column of A is stored
85 in the array AP as follows: if UPLO = 'U', AP(i + (j-1)*j/2) =
86 A(i,j) for 1<=i<=j; if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) =
87 A(i,j) for j<=i<=n. See below for further details.
88 AFP (input or output) COMPLEX*16 array, dimension (N*(N+1)/2)
90 If FACT = 'F', then AFP is an input argument and on entry con‐
91 tains the block diagonal matrix D and the multipliers used to
92 obtain the factor U or L from the factorization A = U*D*U**T or
93 A = L*D*L**T as computed by ZSPTRF, stored as a packed triangu‐
94 lar matrix in the same storage format as A.
95 If FACT = 'N', then AFP is an output argument and on exit con‐
97 tains the block diagonal matrix D and the multipliers used to
98 obtain the factor U or L from the factorization A = U*D*U**T or
99 A = L*D*L**T as computed by ZSPTRF, stored as a packed triangu‐
100 lar matrix in the same storage format as A.
101 IPIV (input or output) INTEGER array, dimension (N)
103 If FACT = 'F', then IPIV is an input argument and on entry con‐
104 tains details of the interchanges and the block structure of D,
105 as determined by ZSPTRF. If IPIV(k) > 0, then rows and columns
106 k and IPIV(k) were interchanged and D(k,k) is a 1-by-1 diagonal
107 block. If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows
108 and columns k-1 and -IPIV(k) were interchanged and
109 D(k-1:k,k-1:k) is a 2-by-2 diagonal block. If UPLO = 'L' and
110 IPIV(k) = IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k)
111 were interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal
112 block.
113 If FACT = 'N', then IPIV is an output argument and on exit con‐
115 tains details of the interchanges and the block structure of D,
116 as determined by ZSPTRF.
117 B (input) COMPLEX*16 array, dimension (LDB,NRHS)
119 The N-by-NRHS right hand side matrix B.
120 LDB (input) INTEGER
122 The leading dimension of the array B. LDB >= max(1,N).
123 X (output) COMPLEX*16 array, dimension (LDX,NRHS)
125 If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
126 LDX (input) INTEGER
128 The leading dimension of the array X. LDX >= max(1,N).
129 RCOND (output) DOUBLE PRECISION
131 The estimate of the reciprocal condition number of the matrix
132 A. If RCOND is less than the machine precision (in particular,
133 if RCOND = 0), the matrix is singular to working precision.
134 This condition is indicated by a return code of INFO > 0.
135 FERR (output) DOUBLE PRECISION array, dimension (NRHS)
137 The estimated forward error bound for each solution vector X(j)
138 (the j-th column of the solution matrix X). If XTRUE is the
139 true solution corresponding to X(j), FERR(j) is an estimated
140 upper bound for the magnitude of the largest element in (X(j) -
141 XTRUE) divided by the magnitude of the largest element in X(j).
142 The estimate is as reliable as the estimate for RCOND, and is
143 almost always a slight overestimate of the true error.
144 BERR (output) DOUBLE PRECISION array, dimension (NRHS)
146 The componentwise relative backward error of each solution vec‐
147 tor X(j) (i.e., the smallest relative change in any element of
148 A or B that makes X(j) an exact solution).
149 WORK (workspace) COMPLEX*16 array, dimension (2*N)
151 RWORK (workspace) DOUBLE PRECISION array, dimension (N)
153 INFO (output) INTEGER
155 = 0: successful exit
156 < 0: if INFO = -i, the i-th argument had an illegal value
157 > 0: if INFO = i, and i is
158 <= N: D(i,i) is exactly zero. The factorization has been com‐
159 pleted but the factor D is exactly singular, so the solution
160 and error bounds could not be computed. RCOND = 0 is returned.
161 = N+1: D is nonsingular, but RCOND is less than machine preci‐
162 sion, meaning that the matrix is singular to working precision.
163 Nevertheless, the solution and error bounds are computed
164 because there are a number of situations where the computed
165 solution can be more accurate than the value of RCOND would
166 suggest.
167
169 The packed storage scheme is illustrated by the following example when
170 N = 4, UPLO = 'U':
171 Two-dimensional storage of the symmetric matrix A:
173 a11 a12 a13 a14
175 a22 a23 a24
176 a33 a34 (aij = aji)
177 a44
178 Packed storage of the upper triangle of A:
180 AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
182 LAPACK driver routine (version 3.N1o)vember 2006 ZSPSVX(1)