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