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