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