1 SUBROUTINE ZSYSVXX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV,
2 $ EQUED, S, B, LDB, X, LDX, RCOND, RPVGRW, BERR,
3 $ N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP,
4 $ NPARAMS, PARAMS, WORK, RWORK, INFO )
5 *
6 * -- LAPACK driver routine (version 3.2.2) --
7 * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
8 * -- Jason Riedy of Univ. of California Berkeley. --
9 * -- June 2010 --
10 *
11 * -- LAPACK is a software package provided by Univ. of Tennessee, --
12 * -- Univ. of California Berkeley and NAG Ltd. --
13 *
14 IMPLICIT NONE
15 * ..
16 * .. Scalar Arguments ..
17 CHARACTER EQUED, FACT, UPLO
18 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
19 $ N_ERR_BNDS
20 DOUBLE PRECISION RCOND, RPVGRW
21 * ..
22 * .. Array Arguments ..
23 INTEGER IPIV( * )
24 COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
25 $ X( LDX, * ), WORK( * )
26 DOUBLE PRECISION S( * ), PARAMS( * ), BERR( * ),
27 $ ERR_BNDS_NORM( NRHS, * ),
28 $ ERR_BNDS_COMP( NRHS, * ), RWORK( * )
29 * ..
30 *
31 * Purpose
32 * =======
33 *
34 * ZSYSVXX uses the diagonal pivoting factorization to compute the
35 * solution to a complex*16 system of linear equations A * X = B, where
36 * A is an N-by-N symmetric matrix and X and B are N-by-NRHS
37 * matrices.
38 *
39 * If requested, both normwise and maximum componentwise error bounds
40 * are returned. ZSYSVXX will return a solution with a tiny
41 * guaranteed error (O(eps) where eps is the working machine
42 * precision) unless the matrix is very ill-conditioned, in which
43 * case a warning is returned. Relevant condition numbers also are
44 * calculated and returned.
45 *
46 * ZSYSVXX accepts user-provided factorizations and equilibration
47 * factors; see the definitions of the FACT and EQUED options.
48 * Solving with refinement and using a factorization from a previous
49 * ZSYSVXX call will also produce a solution with either O(eps)
50 * errors or warnings, but we cannot make that claim for general
51 * user-provided factorizations and equilibration factors if they
52 * differ from what ZSYSVXX would itself produce.
53 *
54 * Description
55 * ===========
56 *
57 * The following steps are performed:
58 *
59 * 1. If FACT = 'E', double precision scaling factors are computed to equilibrate
60 * the system:
61 *
62 * diag(S)*A*diag(S) *inv(diag(S))*X = diag(S)*B
63 *
64 * Whether or not the system will be equilibrated depends on the
65 * scaling of the matrix A, but if equilibration is used, A is
66 * overwritten by diag(S)*A*diag(S) and B by diag(S)*B.
67 *
68 * 2. If FACT = 'N' or 'E', the LU decomposition is used to factor
69 * the matrix A (after equilibration if FACT = 'E') as
70 *
71 * A = U * D * U**T, if UPLO = 'U', or
72 * A = L * D * L**T, if UPLO = 'L',
73 *
74 * where U (or L) is a product of permutation and unit upper (lower)
75 * triangular matrices, and D is symmetric and block diagonal with
76 * 1-by-1 and 2-by-2 diagonal blocks.
77 *
78 * 3. If some D(i,i)=0, so that D is exactly singular, then the
79 * routine returns with INFO = i. Otherwise, the factored form of A
80 * is used to estimate the condition number of the matrix A (see
81 * argument RCOND). If the reciprocal of the condition number is
82 * less than machine precision, the routine still goes on to solve
83 * for X and compute error bounds as described below.
84 *
85 * 4. The system of equations is solved for X using the factored form
86 * of A.
87 *
88 * 5. By default (unless PARAMS(LA_LINRX_ITREF_I) is set to zero),
89 * the routine will use iterative refinement to try to get a small
90 * error and error bounds. Refinement calculates the residual to at
91 * least twice the working precision.
92 *
93 * 6. If equilibration was used, the matrix X is premultiplied by
94 * diag(R) so that it solves the original system before
95 * equilibration.
96 *
97 * Arguments
98 * =========
99 *
100 * Some optional parameters are bundled in the PARAMS array. These
101 * settings determine how refinement is performed, but often the
102 * defaults are acceptable. If the defaults are acceptable, users
103 * can pass NPARAMS = 0 which prevents the source code from accessing
104 * the PARAMS argument.
105 *
106 * FACT (input) CHARACTER*1
107 * Specifies whether or not the factored form of the matrix A is
108 * supplied on entry, and if not, whether the matrix A should be
109 * equilibrated before it is factored.
110 * = 'F': On entry, AF and IPIV contain the factored form of A.
111 * If EQUED is not 'N', the matrix A has been
112 * equilibrated with scaling factors given by S.
113 * A, AF, and IPIV are not modified.
114 * = 'N': The matrix A will be copied to AF and factored.
115 * = 'E': The matrix A will be equilibrated if necessary, then
116 * copied to AF and factored.
117 *
118 * UPLO (input) CHARACTER*1
119 * = 'U': Upper triangle of A is stored;
120 * = 'L': Lower triangle of A is stored.
121 *
122 * N (input) INTEGER
123 * The number of linear equations, i.e., the order of the
124 * matrix A. N >= 0.
125 *
126 * NRHS (input) INTEGER
127 * The number of right hand sides, i.e., the number of columns
128 * of the matrices B and X. NRHS >= 0.
129 *
130 * A (input/output) COMPLEX*16 array, dimension (LDA,N)
131 * The symmetric matrix A. If UPLO = 'U', the leading N-by-N
132 * upper triangular part of A contains the upper triangular
133 * part of the matrix A, and the strictly lower triangular
134 * part of A is not referenced. If UPLO = 'L', the leading
135 * N-by-N lower triangular part of A contains the lower
136 * triangular part of the matrix A, and the strictly upper
137 * triangular part of A is not referenced.
138 *
139 * On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
140 * diag(S)*A*diag(S).
141 *
142 * LDA (input) INTEGER
143 * The leading dimension of the array A. LDA >= max(1,N).
144 *
145 * AF (input or output) COMPLEX*16 array, dimension (LDAF,N)
146 * If FACT = 'F', then AF is an input argument and on entry
147 * contains the block diagonal matrix D and the multipliers
148 * used to obtain the factor U or L from the factorization A =
149 * U*D*U**T or A = L*D*L**T as computed by DSYTRF.
150 *
151 * If FACT = 'N', then AF is an output argument and on exit
152 * returns the block diagonal matrix D and the multipliers
153 * used to obtain the factor U or L from the factorization A =
154 * U*D*U**T or A = L*D*L**T.
155 *
156 * LDAF (input) INTEGER
157 * The leading dimension of the array AF. LDAF >= max(1,N).
158 *
159 * IPIV (input or output) INTEGER array, dimension (N)
160 * If FACT = 'F', then IPIV is an input argument and on entry
161 * contains details of the interchanges and the block
162 * structure of D, as determined by DSYTRF. If IPIV(k) > 0,
163 * then rows and columns k and IPIV(k) were interchanged and
164 * D(k,k) is a 1-by-1 diagonal block. If UPLO = 'U' and
165 * IPIV(k) = IPIV(k-1) < 0, then rows and columns k-1 and
166 * -IPIV(k) were interchanged and D(k-1:k,k-1:k) is a 2-by-2
167 * diagonal block. If UPLO = 'L' and IPIV(k) = IPIV(k+1) < 0,
168 * then rows and columns k+1 and -IPIV(k) were interchanged
169 * and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
170 *
171 * If FACT = 'N', then IPIV is an output argument and on exit
172 * contains details of the interchanges and the block
173 * structure of D, as determined by DSYTRF.
174 *
175 * EQUED (input or output) CHARACTER*1
176 * Specifies the form of equilibration that was done.
177 * = 'N': No equilibration (always true if FACT = 'N').
178 * = 'Y': Both row and column equilibration, i.e., A has been
179 * replaced by diag(S) * A * diag(S).
180 * EQUED is an input argument if FACT = 'F'; otherwise, it is an
181 * output argument.
182 *
183 * S (input or output) DOUBLE PRECISION array, dimension (N)
184 * The scale factors for A. If EQUED = 'Y', A is multiplied on
185 * the left and right by diag(S). S is an input argument if FACT =
186 * 'F'; otherwise, S is an output argument. If FACT = 'F' and EQUED
187 * = 'Y', each element of S must be positive. If S is output, each
188 * element of S is a power of the radix. If S is input, each element
189 * of S should be a power of the radix to ensure a reliable solution
190 * and error estimates. Scaling by powers of the radix does not cause
191 * rounding errors unless the result underflows or overflows.
192 * Rounding errors during scaling lead to refining with a matrix that
193 * is not equivalent to the input matrix, producing error estimates
194 * that may not be reliable.
195 *
196 * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
197 * On entry, the N-by-NRHS right hand side matrix B.
198 * On exit,
199 * if EQUED = 'N', B is not modified;
200 * if EQUED = 'Y', B is overwritten by diag(S)*B;
201 *
202 * LDB (input) INTEGER
203 * The leading dimension of the array B. LDB >= max(1,N).
204 *
205 * X (output) COMPLEX*16 array, dimension (LDX,NRHS)
206 * If INFO = 0, the N-by-NRHS solution matrix X to the original
207 * system of equations. Note that A and B are modified on exit if
208 * EQUED .ne. 'N', and the solution to the equilibrated system is
209 * inv(diag(S))*X.
210 *
211 * LDX (input) INTEGER
212 * The leading dimension of the array X. LDX >= max(1,N).
213 *
214 * RCOND (output) DOUBLE PRECISION
215 * Reciprocal scaled condition number. This is an estimate of the
216 * reciprocal Skeel condition number of the matrix A after
217 * equilibration (if done). If this is less than the machine
218 * precision (in particular, if it is zero), the matrix is singular
219 * to working precision. Note that the error may still be small even
220 * if this number is very small and the matrix appears ill-
221 * conditioned.
222 *
223 * RPVGRW (output) DOUBLE PRECISION
224 * Reciprocal pivot growth. On exit, this contains the reciprocal
225 * pivot growth factor norm(A)/norm(U). The "max absolute element"
226 * norm is used. If this is much less than 1, then the stability of
227 * the LU factorization of the (equilibrated) matrix A could be poor.
228 * This also means that the solution X, estimated condition numbers,
229 * and error bounds could be unreliable. If factorization fails with
230 * 0<INFO<=N, then this contains the reciprocal pivot growth factor
231 * for the leading INFO columns of A.
232 *
233 * BERR (output) DOUBLE PRECISION array, dimension (NRHS)
234 * Componentwise relative backward error. This is the
235 * componentwise relative backward error of each solution vector X(j)
236 * (i.e., the smallest relative change in any element of A or B that
237 * makes X(j) an exact solution).
238 *
239 * N_ERR_BNDS (input) INTEGER
240 * Number of error bounds to return for each right hand side
241 * and each type (normwise or componentwise). See ERR_BNDS_NORM and
242 * ERR_BNDS_COMP below.
243 *
244 * ERR_BNDS_NORM (output) DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
245 * For each right-hand side, this array contains information about
246 * various error bounds and condition numbers corresponding to the
247 * normwise relative error, which is defined as follows:
248 *
249 * Normwise relative error in the ith solution vector:
250 * max_j (abs(XTRUE(j,i) - X(j,i)))
251 * ------------------------------
252 * max_j abs(X(j,i))
253 *
254 * The array is indexed by the type of error information as described
255 * below. There currently are up to three pieces of information
256 * returned.
257 *
258 * The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
259 * right-hand side.
260 *
261 * The second index in ERR_BNDS_NORM(:,err) contains the following
262 * three fields:
263 * err = 1 "Trust/don't trust" boolean. Trust the answer if the
264 * reciprocal condition number is less than the threshold
265 * sqrt(n) * dlamch('Epsilon').
266 *
267 * err = 2 "Guaranteed" error bound: The estimated forward error,
268 * almost certainly within a factor of 10 of the true error
269 * so long as the next entry is greater than the threshold
270 * sqrt(n) * dlamch('Epsilon'). This error bound should only
271 * be trusted if the previous boolean is true.
272 *
273 * err = 3 Reciprocal condition number: Estimated normwise
274 * reciprocal condition number. Compared with the threshold
275 * sqrt(n) * dlamch('Epsilon') to determine if the error
276 * estimate is "guaranteed". These reciprocal condition
277 * numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
278 * appropriately scaled matrix Z.
279 * Let Z = S*A, where S scales each row by a power of the
280 * radix so all absolute row sums of Z are approximately 1.
281 *
282 * See Lapack Working Note 165 for further details and extra
283 * cautions.
284 *
285 * ERR_BNDS_COMP (output) DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
286 * For each right-hand side, this array contains information about
287 * various error bounds and condition numbers corresponding to the
288 * componentwise relative error, which is defined as follows:
289 *
290 * Componentwise relative error in the ith solution vector:
291 * abs(XTRUE(j,i) - X(j,i))
292 * max_j ----------------------
293 * abs(X(j,i))
294 *
295 * The array is indexed by the right-hand side i (on which the
296 * componentwise relative error depends), and the type of error
297 * information as described below. There currently are up to three
298 * pieces of information returned for each right-hand side. If
299 * componentwise accuracy is not requested (PARAMS(3) = 0.0), then
300 * ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most
301 * the first (:,N_ERR_BNDS) entries are returned.
302 *
303 * The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
304 * right-hand side.
305 *
306 * The second index in ERR_BNDS_COMP(:,err) contains the following
307 * three fields:
308 * err = 1 "Trust/don't trust" boolean. Trust the answer if the
309 * reciprocal condition number is less than the threshold
310 * sqrt(n) * dlamch('Epsilon').
311 *
312 * err = 2 "Guaranteed" error bound: The estimated forward error,
313 * almost certainly within a factor of 10 of the true error
314 * so long as the next entry is greater than the threshold
315 * sqrt(n) * dlamch('Epsilon'). This error bound should only
316 * be trusted if the previous boolean is true.
317 *
318 * err = 3 Reciprocal condition number: Estimated componentwise
319 * reciprocal condition number. Compared with the threshold
320 * sqrt(n) * dlamch('Epsilon') to determine if the error
321 * estimate is "guaranteed". These reciprocal condition
322 * numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
323 * appropriately scaled matrix Z.
324 * Let Z = S*(A*diag(x)), where x is the solution for the
325 * current right-hand side and S scales each row of
326 * A*diag(x) by a power of the radix so all absolute row
327 * sums of Z are approximately 1.
328 *
329 * See Lapack Working Note 165 for further details and extra
330 * cautions.
331 *
332 * NPARAMS (input) INTEGER
333 * Specifies the number of parameters set in PARAMS. If .LE. 0, the
334 * PARAMS array is never referenced and default values are used.
335 *
336 * PARAMS (input / output) DOUBLE PRECISION array, dimension NPARAMS
337 * Specifies algorithm parameters. If an entry is .LT. 0.0, then
338 * that entry will be filled with default value used for that
339 * parameter. Only positions up to NPARAMS are accessed; defaults
340 * are used for higher-numbered parameters.
341 *
342 * PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
343 * refinement or not.
344 * Default: 1.0D+0
345 * = 0.0 : No refinement is performed, and no error bounds are
346 * computed.
347 * = 1.0 : Use the extra-precise refinement algorithm.
348 * (other values are reserved for future use)
349 *
350 * PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
351 * computations allowed for refinement.
352 * Default: 10
353 * Aggressive: Set to 100 to permit convergence using approximate
354 * factorizations or factorizations other than LU. If
355 * the factorization uses a technique other than
356 * Gaussian elimination, the guarantees in
357 * err_bnds_norm and err_bnds_comp may no longer be
358 * trustworthy.
359 *
360 * PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
361 * will attempt to find a solution with small componentwise
362 * relative error in the double-precision algorithm. Positive
363 * is true, 0.0 is false.
364 * Default: 1.0 (attempt componentwise convergence)
365 *
366 * WORK (workspace) COMPLEX*16 array, dimension (2*N)
367 *
368 * RWORK (workspace) DOUBLE PRECISION array, dimension (2*N)
369 *
370 * INFO (output) INTEGER
371 * = 0: Successful exit. The solution to every right-hand side is
372 * guaranteed.
373 * < 0: If INFO = -i, the i-th argument had an illegal value
374 * > 0 and <= N: U(INFO,INFO) is exactly zero. The factorization
375 * has been completed, but the factor U is exactly singular, so
376 * the solution and error bounds could not be computed. RCOND = 0
377 * is returned.
378 * = N+J: The solution corresponding to the Jth right-hand side is
379 * not guaranteed. The solutions corresponding to other right-
380 * hand sides K with K > J may not be guaranteed as well, but
381 * only the first such right-hand side is reported. If a small
382 * componentwise error is not requested (PARAMS(3) = 0.0) then
383 * the Jth right-hand side is the first with a normwise error
384 * bound that is not guaranteed (the smallest J such
385 * that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)
386 * the Jth right-hand side is the first with either a normwise or
387 * componentwise error bound that is not guaranteed (the smallest
388 * J such that either ERR_BNDS_NORM(J,1) = 0.0 or
389 * ERR_BNDS_COMP(J,1) = 0.0). See the definition of
390 * ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information
391 * about all of the right-hand sides check ERR_BNDS_NORM or
392 * ERR_BNDS_COMP.
393 *
394 * ==================================================================
395 *
396 * .. Parameters ..
397 DOUBLE PRECISION ZERO, ONE
398 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
399 INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I
400 INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I
401 INTEGER CMP_ERR_I, PIV_GROWTH_I
402 PARAMETER ( FINAL_NRM_ERR_I = 1, FINAL_CMP_ERR_I = 2,
403 $ BERR_I = 3 )
404 PARAMETER ( RCOND_I = 4, NRM_RCOND_I = 5, NRM_ERR_I = 6 )
405 PARAMETER ( CMP_RCOND_I = 7, CMP_ERR_I = 8,
406 $ PIV_GROWTH_I = 9 )
407 * ..
408 * .. Local Scalars ..
409 LOGICAL EQUIL, NOFACT, RCEQU
410 INTEGER INFEQU, J
411 DOUBLE PRECISION AMAX, BIGNUM, SMIN, SMAX, SCOND, SMLNUM
412 * ..
413 * .. External Functions ..
414 EXTERNAL LSAME, DLAMCH, ZLA_SYRPVGRW
415 LOGICAL LSAME
416 DOUBLE PRECISION DLAMCH, ZLA_SYRPVGRW
417 * ..
418 * .. External Subroutines ..
419 EXTERNAL ZSYCON, ZSYEQUB, ZSYTRF, ZSYTRS, ZLACPY,
420 $ ZLAQSY, XERBLA, ZLASCL2, ZSYRFSX
421 * ..
422 * .. Intrinsic Functions ..
423 INTRINSIC MAX, MIN
424 * ..
425 * .. Executable Statements ..
426 *
427 INFO = 0
428 NOFACT = LSAME( FACT, 'N' )
429 EQUIL = LSAME( FACT, 'E' )
430 SMLNUM = DLAMCH( 'Safe minimum' )
431 BIGNUM = ONE / SMLNUM
432 IF( NOFACT .OR. EQUIL ) THEN
433 EQUED = 'N'
434 RCEQU = .FALSE.
435 ELSE
436 RCEQU = LSAME( EQUED, 'Y' )
437 ENDIF
438 *
439 * Default is failure. If an input parameter is wrong or
440 * factorization fails, make everything look horrible. Only the
441 * pivot growth is set here, the rest is initialized in ZSYRFSX.
442 *
443 RPVGRW = ZERO
444 *
445 * Test the input parameters. PARAMS is not tested until ZSYRFSX.
446 *
447 IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.
448 $ LSAME( FACT, 'F' ) ) THEN
449 INFO = -1
450 ELSE IF( .NOT.LSAME(UPLO, 'U') .AND.
451 $ .NOT.LSAME(UPLO, 'L') ) THEN
452 INFO = -2
453 ELSE IF( N.LT.0 ) THEN
454 INFO = -3
455 ELSE IF( NRHS.LT.0 ) THEN
456 INFO = -4
457 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
458 INFO = -6
459 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
460 INFO = -8
461 ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT.
462 $ ( RCEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN
463 INFO = -9
464 ELSE
465 IF ( RCEQU ) THEN
466 SMIN = BIGNUM
467 SMAX = ZERO
468 DO 10 J = 1, N
469 SMIN = MIN( SMIN, S( J ) )
470 SMAX = MAX( SMAX, S( J ) )
471 10 CONTINUE
472 IF( SMIN.LE.ZERO ) THEN
473 INFO = -10
474 ELSE IF( N.GT.0 ) THEN
475 SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM )
476 ELSE
477 SCOND = ONE
478 END IF
479 END IF
480 IF( INFO.EQ.0 ) THEN
481 IF( LDB.LT.MAX( 1, N ) ) THEN
482 INFO = -12
483 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
484 INFO = -14
485 END IF
486 END IF
487 END IF
488 *
489 IF( INFO.NE.0 ) THEN
490 CALL XERBLA( 'ZSYSVXX', -INFO )
491 RETURN
492 END IF
493 *
494 IF( EQUIL ) THEN
495 *
496 * Compute row and column scalings to equilibrate the matrix A.
497 *
498 CALL ZSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFEQU )
499 IF( INFEQU.EQ.0 ) THEN
500 *
501 * Equilibrate the matrix.
502 *
503 CALL ZLAQSY( UPLO, N, A, LDA, S, SCOND, AMAX, EQUED )
504 RCEQU = LSAME( EQUED, 'Y' )
505 END IF
506
507 END IF
508 *
509 * Scale the right hand-side.
510 *
511 IF( RCEQU ) CALL ZLASCL2( N, NRHS, S, B, LDB )
512 *
513 IF( NOFACT .OR. EQUIL ) THEN
514 *
515 * Compute the LDL^T or UDU^T factorization of A.
516 *
517 CALL ZLACPY( UPLO, N, N, A, LDA, AF, LDAF )
518 CALL ZSYTRF( UPLO, N, AF, LDAF, IPIV, WORK, 5*MAX(1,N), INFO )
519 *
520 * Return if INFO is non-zero.
521 *
522 IF( INFO.GT.0 ) THEN
523 *
524 * Pivot in column INFO is exactly 0
525 * Compute the reciprocal pivot growth factor of the
526 * leading rank-deficient INFO columns of A.
527 *
528 IF ( N.GT.0 )
529 $ RPVGRW = ZLA_SYRPVGRW( UPLO, N, INFO, A, LDA, AF,
530 $ LDAF, IPIV, RWORK )
531 RETURN
532 END IF
533 END IF
534 *
535 * Compute the reciprocal pivot growth factor RPVGRW.
536 *
537 IF ( N.GT.0 )
538 $ RPVGRW = ZLA_SYRPVGRW( UPLO, N, INFO, A, LDA, AF, LDAF,
539 $ IPIV, RWORK )
540 *
541 * Compute the solution matrix X.
542 *
543 CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
544 CALL ZSYTRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO )
545 *
546 * Use iterative refinement to improve the computed solution and
547 * compute error bounds and backward error estimates for it.
548 *
549 CALL ZSYRFSX( UPLO, EQUED, N, NRHS, A, LDA, AF, LDAF, IPIV,
550 $ S, B, LDB, X, LDX, RCOND, BERR, N_ERR_BNDS, ERR_BNDS_NORM,
551 $ ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, RWORK, INFO )
552 *
553 * Scale solutions.
554 *
555 IF ( RCEQU ) THEN
556 CALL ZLASCL2 (N, NRHS, S, X, LDX )
557 END IF
558 *
559 RETURN
560 *
561 * End of ZSYSVXX
562 *
563 END
2 $ EQUED, S, B, LDB, X, LDX, RCOND, RPVGRW, BERR,
3 $ N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP,
4 $ NPARAMS, PARAMS, WORK, RWORK, INFO )
5 *
6 * -- LAPACK driver routine (version 3.2.2) --
7 * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
8 * -- Jason Riedy of Univ. of California Berkeley. --
9 * -- June 2010 --
10 *
11 * -- LAPACK is a software package provided by Univ. of Tennessee, --
12 * -- Univ. of California Berkeley and NAG Ltd. --
13 *
14 IMPLICIT NONE
15 * ..
16 * .. Scalar Arguments ..
17 CHARACTER EQUED, FACT, UPLO
18 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
19 $ N_ERR_BNDS
20 DOUBLE PRECISION RCOND, RPVGRW
21 * ..
22 * .. Array Arguments ..
23 INTEGER IPIV( * )
24 COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
25 $ X( LDX, * ), WORK( * )
26 DOUBLE PRECISION S( * ), PARAMS( * ), BERR( * ),
27 $ ERR_BNDS_NORM( NRHS, * ),
28 $ ERR_BNDS_COMP( NRHS, * ), RWORK( * )
29 * ..
30 *
31 * Purpose
32 * =======
33 *
34 * ZSYSVXX uses the diagonal pivoting factorization to compute the
35 * solution to a complex*16 system of linear equations A * X = B, where
36 * A is an N-by-N symmetric matrix and X and B are N-by-NRHS
37 * matrices.
38 *
39 * If requested, both normwise and maximum componentwise error bounds
40 * are returned. ZSYSVXX will return a solution with a tiny
41 * guaranteed error (O(eps) where eps is the working machine
42 * precision) unless the matrix is very ill-conditioned, in which
43 * case a warning is returned. Relevant condition numbers also are
44 * calculated and returned.
45 *
46 * ZSYSVXX accepts user-provided factorizations and equilibration
47 * factors; see the definitions of the FACT and EQUED options.
48 * Solving with refinement and using a factorization from a previous
49 * ZSYSVXX call will also produce a solution with either O(eps)
50 * errors or warnings, but we cannot make that claim for general
51 * user-provided factorizations and equilibration factors if they
52 * differ from what ZSYSVXX would itself produce.
53 *
54 * Description
55 * ===========
56 *
57 * The following steps are performed:
58 *
59 * 1. If FACT = 'E', double precision scaling factors are computed to equilibrate
60 * the system:
61 *
62 * diag(S)*A*diag(S) *inv(diag(S))*X = diag(S)*B
63 *
64 * Whether or not the system will be equilibrated depends on the
65 * scaling of the matrix A, but if equilibration is used, A is
66 * overwritten by diag(S)*A*diag(S) and B by diag(S)*B.
67 *
68 * 2. If FACT = 'N' or 'E', the LU decomposition is used to factor
69 * the matrix A (after equilibration if FACT = 'E') as
70 *
71 * A = U * D * U**T, if UPLO = 'U', or
72 * A = L * D * L**T, if UPLO = 'L',
73 *
74 * where U (or L) is a product of permutation and unit upper (lower)
75 * triangular matrices, and D is symmetric and block diagonal with
76 * 1-by-1 and 2-by-2 diagonal blocks.
77 *
78 * 3. If some D(i,i)=0, so that D is exactly singular, then the
79 * routine returns with INFO = i. Otherwise, the factored form of A
80 * is used to estimate the condition number of the matrix A (see
81 * argument RCOND). If the reciprocal of the condition number is
82 * less than machine precision, the routine still goes on to solve
83 * for X and compute error bounds as described below.
84 *
85 * 4. The system of equations is solved for X using the factored form
86 * of A.
87 *
88 * 5. By default (unless PARAMS(LA_LINRX_ITREF_I) is set to zero),
89 * the routine will use iterative refinement to try to get a small
90 * error and error bounds. Refinement calculates the residual to at
91 * least twice the working precision.
92 *
93 * 6. If equilibration was used, the matrix X is premultiplied by
94 * diag(R) so that it solves the original system before
95 * equilibration.
96 *
97 * Arguments
98 * =========
99 *
100 * Some optional parameters are bundled in the PARAMS array. These
101 * settings determine how refinement is performed, but often the
102 * defaults are acceptable. If the defaults are acceptable, users
103 * can pass NPARAMS = 0 which prevents the source code from accessing
104 * the PARAMS argument.
105 *
106 * FACT (input) CHARACTER*1
107 * Specifies whether or not the factored form of the matrix A is
108 * supplied on entry, and if not, whether the matrix A should be
109 * equilibrated before it is factored.
110 * = 'F': On entry, AF and IPIV contain the factored form of A.
111 * If EQUED is not 'N', the matrix A has been
112 * equilibrated with scaling factors given by S.
113 * A, AF, and IPIV are not modified.
114 * = 'N': The matrix A will be copied to AF and factored.
115 * = 'E': The matrix A will be equilibrated if necessary, then
116 * copied to AF and factored.
117 *
118 * UPLO (input) CHARACTER*1
119 * = 'U': Upper triangle of A is stored;
120 * = 'L': Lower triangle of A is stored.
121 *
122 * N (input) INTEGER
123 * The number of linear equations, i.e., the order of the
124 * matrix A. N >= 0.
125 *
126 * NRHS (input) INTEGER
127 * The number of right hand sides, i.e., the number of columns
128 * of the matrices B and X. NRHS >= 0.
129 *
130 * A (input/output) COMPLEX*16 array, dimension (LDA,N)
131 * The symmetric matrix A. If UPLO = 'U', the leading N-by-N
132 * upper triangular part of A contains the upper triangular
133 * part of the matrix A, and the strictly lower triangular
134 * part of A is not referenced. If UPLO = 'L', the leading
135 * N-by-N lower triangular part of A contains the lower
136 * triangular part of the matrix A, and the strictly upper
137 * triangular part of A is not referenced.
138 *
139 * On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
140 * diag(S)*A*diag(S).
141 *
142 * LDA (input) INTEGER
143 * The leading dimension of the array A. LDA >= max(1,N).
144 *
145 * AF (input or output) COMPLEX*16 array, dimension (LDAF,N)
146 * If FACT = 'F', then AF is an input argument and on entry
147 * contains the block diagonal matrix D and the multipliers
148 * used to obtain the factor U or L from the factorization A =
149 * U*D*U**T or A = L*D*L**T as computed by DSYTRF.
150 *
151 * If FACT = 'N', then AF is an output argument and on exit
152 * returns the block diagonal matrix D and the multipliers
153 * used to obtain the factor U or L from the factorization A =
154 * U*D*U**T or A = L*D*L**T.
155 *
156 * LDAF (input) INTEGER
157 * The leading dimension of the array AF. LDAF >= max(1,N).
158 *
159 * IPIV (input or output) INTEGER array, dimension (N)
160 * If FACT = 'F', then IPIV is an input argument and on entry
161 * contains details of the interchanges and the block
162 * structure of D, as determined by DSYTRF. If IPIV(k) > 0,
163 * then rows and columns k and IPIV(k) were interchanged and
164 * D(k,k) is a 1-by-1 diagonal block. If UPLO = 'U' and
165 * IPIV(k) = IPIV(k-1) < 0, then rows and columns k-1 and
166 * -IPIV(k) were interchanged and D(k-1:k,k-1:k) is a 2-by-2
167 * diagonal block. If UPLO = 'L' and IPIV(k) = IPIV(k+1) < 0,
168 * then rows and columns k+1 and -IPIV(k) were interchanged
169 * and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
170 *
171 * If FACT = 'N', then IPIV is an output argument and on exit
172 * contains details of the interchanges and the block
173 * structure of D, as determined by DSYTRF.
174 *
175 * EQUED (input or output) CHARACTER*1
176 * Specifies the form of equilibration that was done.
177 * = 'N': No equilibration (always true if FACT = 'N').
178 * = 'Y': Both row and column equilibration, i.e., A has been
179 * replaced by diag(S) * A * diag(S).
180 * EQUED is an input argument if FACT = 'F'; otherwise, it is an
181 * output argument.
182 *
183 * S (input or output) DOUBLE PRECISION array, dimension (N)
184 * The scale factors for A. If EQUED = 'Y', A is multiplied on
185 * the left and right by diag(S). S is an input argument if FACT =
186 * 'F'; otherwise, S is an output argument. If FACT = 'F' and EQUED
187 * = 'Y', each element of S must be positive. If S is output, each
188 * element of S is a power of the radix. If S is input, each element
189 * of S should be a power of the radix to ensure a reliable solution
190 * and error estimates. Scaling by powers of the radix does not cause
191 * rounding errors unless the result underflows or overflows.
192 * Rounding errors during scaling lead to refining with a matrix that
193 * is not equivalent to the input matrix, producing error estimates
194 * that may not be reliable.
195 *
196 * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
197 * On entry, the N-by-NRHS right hand side matrix B.
198 * On exit,
199 * if EQUED = 'N', B is not modified;
200 * if EQUED = 'Y', B is overwritten by diag(S)*B;
201 *
202 * LDB (input) INTEGER
203 * The leading dimension of the array B. LDB >= max(1,N).
204 *
205 * X (output) COMPLEX*16 array, dimension (LDX,NRHS)
206 * If INFO = 0, the N-by-NRHS solution matrix X to the original
207 * system of equations. Note that A and B are modified on exit if
208 * EQUED .ne. 'N', and the solution to the equilibrated system is
209 * inv(diag(S))*X.
210 *
211 * LDX (input) INTEGER
212 * The leading dimension of the array X. LDX >= max(1,N).
213 *
214 * RCOND (output) DOUBLE PRECISION
215 * Reciprocal scaled condition number. This is an estimate of the
216 * reciprocal Skeel condition number of the matrix A after
217 * equilibration (if done). If this is less than the machine
218 * precision (in particular, if it is zero), the matrix is singular
219 * to working precision. Note that the error may still be small even
220 * if this number is very small and the matrix appears ill-
221 * conditioned.
222 *
223 * RPVGRW (output) DOUBLE PRECISION
224 * Reciprocal pivot growth. On exit, this contains the reciprocal
225 * pivot growth factor norm(A)/norm(U). The "max absolute element"
226 * norm is used. If this is much less than 1, then the stability of
227 * the LU factorization of the (equilibrated) matrix A could be poor.
228 * This also means that the solution X, estimated condition numbers,
229 * and error bounds could be unreliable. If factorization fails with
230 * 0<INFO<=N, then this contains the reciprocal pivot growth factor
231 * for the leading INFO columns of A.
232 *
233 * BERR (output) DOUBLE PRECISION array, dimension (NRHS)
234 * Componentwise relative backward error. This is the
235 * componentwise relative backward error of each solution vector X(j)
236 * (i.e., the smallest relative change in any element of A or B that
237 * makes X(j) an exact solution).
238 *
239 * N_ERR_BNDS (input) INTEGER
240 * Number of error bounds to return for each right hand side
241 * and each type (normwise or componentwise). See ERR_BNDS_NORM and
242 * ERR_BNDS_COMP below.
243 *
244 * ERR_BNDS_NORM (output) DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
245 * For each right-hand side, this array contains information about
246 * various error bounds and condition numbers corresponding to the
247 * normwise relative error, which is defined as follows:
248 *
249 * Normwise relative error in the ith solution vector:
250 * max_j (abs(XTRUE(j,i) - X(j,i)))
251 * ------------------------------
252 * max_j abs(X(j,i))
253 *
254 * The array is indexed by the type of error information as described
255 * below. There currently are up to three pieces of information
256 * returned.
257 *
258 * The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
259 * right-hand side.
260 *
261 * The second index in ERR_BNDS_NORM(:,err) contains the following
262 * three fields:
263 * err = 1 "Trust/don't trust" boolean. Trust the answer if the
264 * reciprocal condition number is less than the threshold
265 * sqrt(n) * dlamch('Epsilon').
266 *
267 * err = 2 "Guaranteed" error bound: The estimated forward error,
268 * almost certainly within a factor of 10 of the true error
269 * so long as the next entry is greater than the threshold
270 * sqrt(n) * dlamch('Epsilon'). This error bound should only
271 * be trusted if the previous boolean is true.
272 *
273 * err = 3 Reciprocal condition number: Estimated normwise
274 * reciprocal condition number. Compared with the threshold
275 * sqrt(n) * dlamch('Epsilon') to determine if the error
276 * estimate is "guaranteed". These reciprocal condition
277 * numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
278 * appropriately scaled matrix Z.
279 * Let Z = S*A, where S scales each row by a power of the
280 * radix so all absolute row sums of Z are approximately 1.
281 *
282 * See Lapack Working Note 165 for further details and extra
283 * cautions.
284 *
285 * ERR_BNDS_COMP (output) DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
286 * For each right-hand side, this array contains information about
287 * various error bounds and condition numbers corresponding to the
288 * componentwise relative error, which is defined as follows:
289 *
290 * Componentwise relative error in the ith solution vector:
291 * abs(XTRUE(j,i) - X(j,i))
292 * max_j ----------------------
293 * abs(X(j,i))
294 *
295 * The array is indexed by the right-hand side i (on which the
296 * componentwise relative error depends), and the type of error
297 * information as described below. There currently are up to three
298 * pieces of information returned for each right-hand side. If
299 * componentwise accuracy is not requested (PARAMS(3) = 0.0), then
300 * ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most
301 * the first (:,N_ERR_BNDS) entries are returned.
302 *
303 * The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
304 * right-hand side.
305 *
306 * The second index in ERR_BNDS_COMP(:,err) contains the following
307 * three fields:
308 * err = 1 "Trust/don't trust" boolean. Trust the answer if the
309 * reciprocal condition number is less than the threshold
310 * sqrt(n) * dlamch('Epsilon').
311 *
312 * err = 2 "Guaranteed" error bound: The estimated forward error,
313 * almost certainly within a factor of 10 of the true error
314 * so long as the next entry is greater than the threshold
315 * sqrt(n) * dlamch('Epsilon'). This error bound should only
316 * be trusted if the previous boolean is true.
317 *
318 * err = 3 Reciprocal condition number: Estimated componentwise
319 * reciprocal condition number. Compared with the threshold
320 * sqrt(n) * dlamch('Epsilon') to determine if the error
321 * estimate is "guaranteed". These reciprocal condition
322 * numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
323 * appropriately scaled matrix Z.
324 * Let Z = S*(A*diag(x)), where x is the solution for the
325 * current right-hand side and S scales each row of
326 * A*diag(x) by a power of the radix so all absolute row
327 * sums of Z are approximately 1.
328 *
329 * See Lapack Working Note 165 for further details and extra
330 * cautions.
331 *
332 * NPARAMS (input) INTEGER
333 * Specifies the number of parameters set in PARAMS. If .LE. 0, the
334 * PARAMS array is never referenced and default values are used.
335 *
336 * PARAMS (input / output) DOUBLE PRECISION array, dimension NPARAMS
337 * Specifies algorithm parameters. If an entry is .LT. 0.0, then
338 * that entry will be filled with default value used for that
339 * parameter. Only positions up to NPARAMS are accessed; defaults
340 * are used for higher-numbered parameters.
341 *
342 * PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
343 * refinement or not.
344 * Default: 1.0D+0
345 * = 0.0 : No refinement is performed, and no error bounds are
346 * computed.
347 * = 1.0 : Use the extra-precise refinement algorithm.
348 * (other values are reserved for future use)
349 *
350 * PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
351 * computations allowed for refinement.
352 * Default: 10
353 * Aggressive: Set to 100 to permit convergence using approximate
354 * factorizations or factorizations other than LU. If
355 * the factorization uses a technique other than
356 * Gaussian elimination, the guarantees in
357 * err_bnds_norm and err_bnds_comp may no longer be
358 * trustworthy.
359 *
360 * PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
361 * will attempt to find a solution with small componentwise
362 * relative error in the double-precision algorithm. Positive
363 * is true, 0.0 is false.
364 * Default: 1.0 (attempt componentwise convergence)
365 *
366 * WORK (workspace) COMPLEX*16 array, dimension (2*N)
367 *
368 * RWORK (workspace) DOUBLE PRECISION array, dimension (2*N)
369 *
370 * INFO (output) INTEGER
371 * = 0: Successful exit. The solution to every right-hand side is
372 * guaranteed.
373 * < 0: If INFO = -i, the i-th argument had an illegal value
374 * > 0 and <= N: U(INFO,INFO) is exactly zero. The factorization
375 * has been completed, but the factor U is exactly singular, so
376 * the solution and error bounds could not be computed. RCOND = 0
377 * is returned.
378 * = N+J: The solution corresponding to the Jth right-hand side is
379 * not guaranteed. The solutions corresponding to other right-
380 * hand sides K with K > J may not be guaranteed as well, but
381 * only the first such right-hand side is reported. If a small
382 * componentwise error is not requested (PARAMS(3) = 0.0) then
383 * the Jth right-hand side is the first with a normwise error
384 * bound that is not guaranteed (the smallest J such
385 * that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)
386 * the Jth right-hand side is the first with either a normwise or
387 * componentwise error bound that is not guaranteed (the smallest
388 * J such that either ERR_BNDS_NORM(J,1) = 0.0 or
389 * ERR_BNDS_COMP(J,1) = 0.0). See the definition of
390 * ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information
391 * about all of the right-hand sides check ERR_BNDS_NORM or
392 * ERR_BNDS_COMP.
393 *
394 * ==================================================================
395 *
396 * .. Parameters ..
397 DOUBLE PRECISION ZERO, ONE
398 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
399 INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I
400 INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I
401 INTEGER CMP_ERR_I, PIV_GROWTH_I
402 PARAMETER ( FINAL_NRM_ERR_I = 1, FINAL_CMP_ERR_I = 2,
403 $ BERR_I = 3 )
404 PARAMETER ( RCOND_I = 4, NRM_RCOND_I = 5, NRM_ERR_I = 6 )
405 PARAMETER ( CMP_RCOND_I = 7, CMP_ERR_I = 8,
406 $ PIV_GROWTH_I = 9 )
407 * ..
408 * .. Local Scalars ..
409 LOGICAL EQUIL, NOFACT, RCEQU
410 INTEGER INFEQU, J
411 DOUBLE PRECISION AMAX, BIGNUM, SMIN, SMAX, SCOND, SMLNUM
412 * ..
413 * .. External Functions ..
414 EXTERNAL LSAME, DLAMCH, ZLA_SYRPVGRW
415 LOGICAL LSAME
416 DOUBLE PRECISION DLAMCH, ZLA_SYRPVGRW
417 * ..
418 * .. External Subroutines ..
419 EXTERNAL ZSYCON, ZSYEQUB, ZSYTRF, ZSYTRS, ZLACPY,
420 $ ZLAQSY, XERBLA, ZLASCL2, ZSYRFSX
421 * ..
422 * .. Intrinsic Functions ..
423 INTRINSIC MAX, MIN
424 * ..
425 * .. Executable Statements ..
426 *
427 INFO = 0
428 NOFACT = LSAME( FACT, 'N' )
429 EQUIL = LSAME( FACT, 'E' )
430 SMLNUM = DLAMCH( 'Safe minimum' )
431 BIGNUM = ONE / SMLNUM
432 IF( NOFACT .OR. EQUIL ) THEN
433 EQUED = 'N'
434 RCEQU = .FALSE.
435 ELSE
436 RCEQU = LSAME( EQUED, 'Y' )
437 ENDIF
438 *
439 * Default is failure. If an input parameter is wrong or
440 * factorization fails, make everything look horrible. Only the
441 * pivot growth is set here, the rest is initialized in ZSYRFSX.
442 *
443 RPVGRW = ZERO
444 *
445 * Test the input parameters. PARAMS is not tested until ZSYRFSX.
446 *
447 IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.
448 $ LSAME( FACT, 'F' ) ) THEN
449 INFO = -1
450 ELSE IF( .NOT.LSAME(UPLO, 'U') .AND.
451 $ .NOT.LSAME(UPLO, 'L') ) THEN
452 INFO = -2
453 ELSE IF( N.LT.0 ) THEN
454 INFO = -3
455 ELSE IF( NRHS.LT.0 ) THEN
456 INFO = -4
457 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
458 INFO = -6
459 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
460 INFO = -8
461 ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT.
462 $ ( RCEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN
463 INFO = -9
464 ELSE
465 IF ( RCEQU ) THEN
466 SMIN = BIGNUM
467 SMAX = ZERO
468 DO 10 J = 1, N
469 SMIN = MIN( SMIN, S( J ) )
470 SMAX = MAX( SMAX, S( J ) )
471 10 CONTINUE
472 IF( SMIN.LE.ZERO ) THEN
473 INFO = -10
474 ELSE IF( N.GT.0 ) THEN
475 SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM )
476 ELSE
477 SCOND = ONE
478 END IF
479 END IF
480 IF( INFO.EQ.0 ) THEN
481 IF( LDB.LT.MAX( 1, N ) ) THEN
482 INFO = -12
483 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
484 INFO = -14
485 END IF
486 END IF
487 END IF
488 *
489 IF( INFO.NE.0 ) THEN
490 CALL XERBLA( 'ZSYSVXX', -INFO )
491 RETURN
492 END IF
493 *
494 IF( EQUIL ) THEN
495 *
496 * Compute row and column scalings to equilibrate the matrix A.
497 *
498 CALL ZSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFEQU )
499 IF( INFEQU.EQ.0 ) THEN
500 *
501 * Equilibrate the matrix.
502 *
503 CALL ZLAQSY( UPLO, N, A, LDA, S, SCOND, AMAX, EQUED )
504 RCEQU = LSAME( EQUED, 'Y' )
505 END IF
506
507 END IF
508 *
509 * Scale the right hand-side.
510 *
511 IF( RCEQU ) CALL ZLASCL2( N, NRHS, S, B, LDB )
512 *
513 IF( NOFACT .OR. EQUIL ) THEN
514 *
515 * Compute the LDL^T or UDU^T factorization of A.
516 *
517 CALL ZLACPY( UPLO, N, N, A, LDA, AF, LDAF )
518 CALL ZSYTRF( UPLO, N, AF, LDAF, IPIV, WORK, 5*MAX(1,N), INFO )
519 *
520 * Return if INFO is non-zero.
521 *
522 IF( INFO.GT.0 ) THEN
523 *
524 * Pivot in column INFO is exactly 0
525 * Compute the reciprocal pivot growth factor of the
526 * leading rank-deficient INFO columns of A.
527 *
528 IF ( N.GT.0 )
529 $ RPVGRW = ZLA_SYRPVGRW( UPLO, N, INFO, A, LDA, AF,
530 $ LDAF, IPIV, RWORK )
531 RETURN
532 END IF
533 END IF
534 *
535 * Compute the reciprocal pivot growth factor RPVGRW.
536 *
537 IF ( N.GT.0 )
538 $ RPVGRW = ZLA_SYRPVGRW( UPLO, N, INFO, A, LDA, AF, LDAF,
539 $ IPIV, RWORK )
540 *
541 * Compute the solution matrix X.
542 *
543 CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
544 CALL ZSYTRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO )
545 *
546 * Use iterative refinement to improve the computed solution and
547 * compute error bounds and backward error estimates for it.
548 *
549 CALL ZSYRFSX( UPLO, EQUED, N, NRHS, A, LDA, AF, LDAF, IPIV,
550 $ S, B, LDB, X, LDX, RCOND, BERR, N_ERR_BNDS, ERR_BNDS_NORM,
551 $ ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, RWORK, INFO )
552 *
553 * Scale solutions.
554 *
555 IF ( RCEQU ) THEN
556 CALL ZLASCL2 (N, NRHS, S, X, LDX )
557 END IF
558 *
559 RETURN
560 *
561 * End of ZSYSVXX
562 *
563 END