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