1ZSYRFSX(1) LAPACK routine (version 3.2) ZSYRFSX(1)
2
6 ZSYRFSX - ZSYRFSX improve the computed solution to a system of linear
7 equations when the coefficient matrix is symmetric indefinite, and
8 provides error bounds and backward error estimates for the solution
9
11 Subroutine ZSYRFSX( UPLO, EQUED, N, NRHS, A, LDA, AF, LDAF, IPIV, S, B,
12 LDB, X, LDX, RCOND, BERR, N_ERR_BNDS,
13 ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS,
14 WORK, RWORK, INFO )
15 IMPLICIT NONE
17 CHARACTER UPLO, EQUED
19 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
21 N_ERR_BNDS
22 DOUBLE PRECISION RCOND
24 INTEGER IPIV( * )
26 COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ), X( LDX, *
28 ), WORK( * )
29 DOUBLE PRECISION S( * ), PARAMS( * ), BERR( * ), RWORK( *
31 ), ERR_BNDS_NORM( NRHS, * ), ERR_BNDS_COMP( NRHS, *
32 )
33
35 ZSYRFSX improves the computed solution to a system of linear
36 equations when the coefficient matrix is symmetric indefinite, and
37 provides error bounds and backward error estimates for the
38 solution. In addition to normwise error bound, the code provides
39 maximum componentwise error bound if possible. See comments for
40 ERR_BNDS_N and ERR_BNDS_C for details of the error bounds.
41 The original system of linear equations may have been equilibrated
42 before calling this routine, as described by arguments EQUED and S
43 below. In this case, the solution and error bounds returned are
44 for the original unequilibrated system.
45
47 Some optional parameters are bundled in the PARAMS array. These set‐
48 tings determine how refinement is performed, but often the defaults are
49 acceptable. If the defaults are acceptable, users can pass NPARAMS = 0
50 which prevents the source code from accessing the PARAMS argument.
51 UPLO (input) CHARACTER*1
53 = 'U': Upper triangle of A is stored;
54 = 'L': Lower triangle of A is stored.
55 EQUED (input) CHARACTER*1
57 Specifies the form of equilibration that was done to A before
58 calling this routine. This is needed to compute the solution
59 and error bounds correctly. = 'N': No equilibration
60 = 'Y': Both row and column equilibration, i.e., A has been
61 replaced by diag(S) * A * diag(S). The right hand side B has
62 been changed accordingly.
63 N (input) INTEGER
65 The order of the matrix A. N >= 0.
66 NRHS (input) INTEGER
68 The number of right hand sides, i.e., the number of columns of
69 the matrices B and X. NRHS >= 0.
70 A (input) COMPLEX*16 array, dimension (LDA,N)
72 The symmetric matrix A. If UPLO = 'U', the leading N-by-N
73 upper triangular part of A contains the upper triangular part
74 of the matrix A, and the strictly lower triangular part of A is
75 not referenced. If UPLO = 'L', the leading N-by-N lower trian‐
76 gular part of A contains the lower triangular part of the
77 matrix A, and the strictly upper triangular part of A is not
78 referenced.
79 LDA (input) INTEGER
81 The leading dimension of the array A. LDA >= max(1,N).
82 AF (input) COMPLEX*16 array, dimension (LDAF,N)
84 The factored form of the matrix A. AF contains the block diag‐
85 onal matrix D and the multipliers used to obtain the factor U
86 or L from the factorization A = U*D*U**T or A = L*D*L**T as
87 computed by DSYTRF.
88 LDAF (input) INTEGER
90 The leading dimension of the array AF. LDAF >= max(1,N).
91 IPIV (input) INTEGER array, dimension (N)
93 Details of the interchanges and the block structure of D as
94 determined by DSYTRF.
95 S (input or output) DOUBLE PRECISION array, dimension (N)
97 The scale factors for A. If EQUED = 'Y', A is multiplied on
98 the left and right by diag(S). S is an input argument if FACT
99 = 'F'; otherwise, S is an output argument. If FACT = 'F' and
100 EQUED = 'Y', each element of S must be positive. If S is out‐
101 put, each element of S is a power of the radix. If S is input,
102 each element of S should be a power of the radix to ensure a
103 reliable solution and error estimates. Scaling by powers of the
104 radix does not cause rounding errors unless the result under‐
105 flows or overflows. Rounding errors during scaling lead to
106 refining with a matrix that is not equivalent to the input
107 matrix, producing error estimates that may not be reliable.
108 B (input) COMPLEX*16 array, dimension (LDB,NRHS)
110 The right hand side matrix B.
111 LDB (input) INTEGER
113 The leading dimension of the array B. LDB >= max(1,N).
114 X (input/output) COMPLEX*16 array, dimension (LDX,NRHS)
116 On entry, the solution matrix X, as computed by DGETRS. On
117 exit, the improved solution matrix X.
118 LDX (input) INTEGER
120 The leading dimension of the array X. LDX >= max(1,N).
121 RCOND (output) DOUBLE PRECISION
123 Reciprocal scaled condition number. This is an estimate of the
124 reciprocal Skeel condition number of the matrix A after equili‐
125 bration (if done). If this is less than the machine precision
126 (in particular, if it is zero), the matrix is singular to work‐
127 ing precision. Note that the error may still be small even if
128 this number is very small and the matrix appears ill- condi‐
129 tioned.
130 BERR (output) DOUBLE PRECISION array, dimension (NRHS)
132 Componentwise relative backward error. This is the component‐
133 wise relative backward error of each solution vector X(j)
134 (i.e., the smallest relative change in any element of A or B
135 that makes X(j) an exact solution). N_ERR_BNDS (input) INTEGER
136 Number of error bounds to return for each right hand side and
137 each type (normwise or componentwise). See ERR_BNDS_NORM and
138 ERR_BNDS_COMP below.
139 ERR_BNDS_NORM (output) DOUBLE PRECISION array, dimension (NRHS,
141 N_ERR_BNDS)
142 For each right-hand side, this array contains informa‐
143 tion about various error bounds and condition numbers
144 corresponding to the normwise relative error, which is
145 defined as follows: Normwise relative error in the ith
146 solution vector: max_j (abs(XTRUE(j,i) - X(j,i)))
147 ------------------------------ max_j abs(X(j,i)) The
148 array is indexed by the type of error information as
149 described below. There currently are up to three pieces
150 of information returned. The first index in
151 ERR_BNDS_NORM(i,:) corresponds to the ith right-hand
152 side. The second index in ERR_BNDS_NORM(:,err) contains
153 the following three fields: err = 1 "Trust/don't trust"
154 boolean. Trust the answer if the reciprocal condition
155 number is less than the threshold sqrt(n) *
156 dlamch('Epsilon'). err = 2 "Guaranteed" error bound:
157 The estimated forward error, almost certainly within a
158 factor of 10 of the true error so long as the next entry
159 is greater than the threshold sqrt(n) *
160 dlamch('Epsilon'). This error bound should only be
161 trusted if the previous boolean is true. err = 3
162 Reciprocal condition number: Estimated normwise recipro‐
163 cal condition number. Compared with the threshold
164 sqrt(n) * dlamch('Epsilon') to determine if the error
165 estimate is "guaranteed". These reciprocal condition
166 numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for
167 some appropriately scaled matrix Z. Let Z = S*A, where
168 S scales each row by a power of the radix so all abso‐
169 lute row sums of Z are approximately 1. See Lapack
170 Working Note 165 for further details and extra cautions.
171 ERR_BNDS_COMP (output) DOUBLE PRECISION array, dimension (NRHS,
173 N_ERR_BNDS)
174 For each right-hand side, this array contains informa‐
175 tion about various error bounds and condition numbers
176 corresponding to the componentwise relative error, which
177 is defined as follows: Componentwise relative error in
178 the ith solution vector: abs(XTRUE(j,i) - X(j,i)) max_j
179 ---------------------- abs(X(j,i)) The array is indexed
180 by the right-hand side i (on which the componentwise
181 relative error depends), and the type of error informa‐
182 tion as described below. There currently are up to three
183 pieces of information returned for each right-hand side.
184 If componentwise accuracy is not requested (PARAMS(3) =
185 0.0), then ERR_BNDS_COMP is not accessed. If N_ERR_BNDS
186 .LT. 3, then at most the first (:,N_ERR_BNDS) entries
187 are returned. The first index in ERR_BNDS_COMP(i,:)
188 corresponds to the ith right-hand side. The second
189 index in ERR_BNDS_COMP(:,err) contains the following
190 three fields: err = 1 "Trust/don't trust" boolean. Trust
191 the answer if the reciprocal condition number is less
192 than the threshold sqrt(n) * dlamch('Epsilon'). err = 2
193 "Guaranteed" error bound: The estimated forward error,
194 almost certainly within a factor of 10 of the true error
195 so long as the next entry is greater than the threshold
196 sqrt(n) * dlamch('Epsilon'). This error bound should
197 only be trusted if the previous boolean is true. err =
198 3 Reciprocal condition number: Estimated componentwise
199 reciprocal condition number. Compared with the thresh‐
200 old sqrt(n) * dlamch('Epsilon') to determine if the
201 error estimate is "guaranteed". These reciprocal condi‐
202 tion numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf))
203 for some appropriately scaled matrix Z. Let Z =
204 S*(A*diag(x)), where x is the solution for the current
205 right-hand side and S scales each row of A*diag(x) by a
206 power of the radix so all absolute row sums of Z are
207 approximately 1. See Lapack Working Note 165 for fur‐
208 ther details and extra cautions. NPARAMS (input) INTE‐
209 GER Specifies the number of parameters set in PARAMS.
210 If .LE. 0, the PARAMS array is never referenced and
211 default values are used.
212 PARAMS (input / output) DOUBLE PRECISION array, dimension NPARAMS
214 Specifies algorithm parameters. If an entry is .LT. 0.0, then
215 that entry will be filled with default value used for that
216 parameter. Only positions up to NPARAMS are accessed; defaults
217 are used for higher-numbered parameters.
218 PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
219 refinement or not. Default: 1.0D+0
220 = 0.0 : No refinement is performed, and no error bounds are
221 computed. = 1.0 : Use the double-precision refinement algo‐
222 rithm, possibly with doubled-single computations if the compi‐
223 lation environment does not support DOUBLE PRECISION. (other
224 values are reserved for future use) PARAMS(LA_LINRX_ITHRESH_I =
225 2) : Maximum number of residual computations allowed for
226 refinement. Default: 10
227 Aggressive: Set to 100 to permit convergence using approximate
228 factorizations or factorizations other than LU. If the factor‐
229 ization uses a technique other than Gaussian elimination, the
230 guarantees in err_bnds_norm and err_bnds_comp may no longer be
231 trustworthy. PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining
232 if the code will attempt to find a solution with small compo‐
233 nentwise relative error in the double-precision algorithm.
234 Positive is true, 0.0 is false. Default: 1.0 (attempt compo‐
235 nentwise convergence)
236 WORK (workspace) DOUBLE PRECISION array, dimension (4*N)
238 IWORK (workspace) INTEGER array, dimension (N)
240 INFO (output) INTEGER
242 = 0: Successful exit. The solution to every right-hand side is
243 guaranteed. < 0: If INFO = -i, the i-th argument had an ille‐
244 gal value
245 > 0 and <= N: U(INFO,INFO) is exactly zero. The factorization
246 has been completed, but the factor U is exactly singular, so
247 the solution and error bounds could not be computed. RCOND = 0
248 is returned. = N+J: The solution corresponding to the Jth
249 right-hand side is not guaranteed. The solutions corresponding
250 to other right- hand sides K with K > J may not be guaranteed
251 as well, but only the first such right-hand side is reported.
252 If a small componentwise error is not requested (PARAMS(3) =
253 0.0) then the Jth right-hand side is the first with a normwise
254 error bound that is not guaranteed (the smallest J such that
255 ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0) the Jth
256 right-hand side is the first with either a normwise or compo‐
257 nentwise error bound that is not guaranteed (the smallest J
258 such that either ERR_BNDS_NORM(J,1) = 0.0 or ERR_BNDS_COMP(J,1)
259 = 0.0). See the definition of ERR_BNDS_NORM(:,1) and
260 ERR_BNDS_COMP(:,1). To get information about all of the right-
261 hand sides check ERR_BNDS_NORM or ERR_BNDS_COMP.
262 LAPACK routine (version 3.2) November 2008 ZSYRFSX(1)