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