1*> \brief <b> ZPOSVX computes the solution to system of linear equations A * X = B for PO matrices</b> 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZPOSVX + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zposvx.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zposvx.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zposvx.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZPOSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED, 22* S, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, 23* RWORK, INFO ) 24* 25* .. Scalar Arguments .. 26* CHARACTER EQUED, FACT, UPLO 27* INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS 28* DOUBLE PRECISION RCOND 29* .. 30* .. Array Arguments .. 31* DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * ), S( * ) 32* COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 33* $ WORK( * ), X( LDX, * ) 34* .. 35* 36* 37*> \par Purpose: 38* ============= 39*> 40*> \verbatim 41*> 42*> ZPOSVX uses the Cholesky factorization A = U**H*U or A = L*L**H to 43*> compute the solution to a complex system of linear equations 44*> A * X = B, 45*> where A is an N-by-N Hermitian positive definite matrix and X and B 46*> are N-by-NRHS matrices. 47*> 48*> Error bounds on the solution and a condition estimate are also 49*> provided. 50*> \endverbatim 51* 52*> \par Description: 53* ================= 54*> 55*> \verbatim 56*> 57*> The following steps are performed: 58*> 59*> 1. If FACT = 'E', real scaling factors are computed to equilibrate 60*> the system: 61*> diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B 62*> Whether or not the system will be equilibrated depends on the 63*> scaling of the matrix A, but if equilibration is used, A is 64*> overwritten by diag(S)*A*diag(S) and B by diag(S)*B. 65*> 66*> 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to 67*> factor the matrix A (after equilibration if FACT = 'E') as 68*> A = U**H* U, if UPLO = 'U', or 69*> A = L * L**H, if UPLO = 'L', 70*> where U is an upper triangular matrix and L is a lower triangular 71*> matrix. 72*> 73*> 3. If the leading i-by-i principal minor is not positive definite, 74*> then the routine returns with INFO = i. Otherwise, the factored 75*> form of A is used to estimate the condition number of the matrix 76*> A. If the reciprocal of the condition number is less than machine 77*> precision, INFO = N+1 is returned as a warning, but the routine 78*> still goes on to solve for X and compute error bounds as 79*> described below. 80*> 81*> 4. The system of equations is solved for X using the factored form 82*> of A. 83*> 84*> 5. Iterative refinement is applied to improve the computed solution 85*> matrix and calculate error bounds and backward error estimates 86*> for it. 87*> 88*> 6. If equilibration was used, the matrix X is premultiplied by 89*> diag(S) so that it solves the original system before 90*> equilibration. 91*> \endverbatim 92* 93* Arguments: 94* ========== 95* 96*> \param[in] FACT 97*> \verbatim 98*> FACT is CHARACTER*1 99*> Specifies whether or not the factored form of the matrix A is 100*> supplied on entry, and if not, whether the matrix A should be 101*> equilibrated before it is factored. 102*> = 'F': On entry, AF contains the factored form of A. 103*> If EQUED = 'Y', the matrix A has been equilibrated 104*> with scaling factors given by S. A and AF will not 105*> be modified. 106*> = 'N': The matrix A will be copied to AF and factored. 107*> = 'E': The matrix A will be equilibrated if necessary, then 108*> copied to AF and factored. 109*> \endverbatim 110*> 111*> \param[in] UPLO 112*> \verbatim 113*> UPLO is CHARACTER*1 114*> = 'U': Upper triangle of A is stored; 115*> = 'L': Lower triangle of A is stored. 116*> \endverbatim 117*> 118*> \param[in] N 119*> \verbatim 120*> N is INTEGER 121*> The number of linear equations, i.e., the order of the 122*> matrix A. N >= 0. 123*> \endverbatim 124*> 125*> \param[in] NRHS 126*> \verbatim 127*> NRHS is INTEGER 128*> The number of right hand sides, i.e., the number of columns 129*> of the matrices B and X. NRHS >= 0. 130*> \endverbatim 131*> 132*> \param[in,out] A 133*> \verbatim 134*> A is COMPLEX*16 array, dimension (LDA,N) 135*> On entry, the Hermitian matrix A, except if FACT = 'F' and 136*> EQUED = 'Y', then A must contain the equilibrated matrix 137*> diag(S)*A*diag(S). If UPLO = 'U', the leading 138*> N-by-N upper triangular part of A contains the upper 139*> triangular part of the matrix A, and the strictly lower 140*> triangular part of A is not referenced. If UPLO = 'L', the 141*> leading N-by-N lower triangular part of A contains the lower 142*> triangular part of the matrix A, and the strictly upper 143*> triangular part of A is not referenced. A is not modified if 144*> FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit. 145*> 146*> On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by 147*> diag(S)*A*diag(S). 148*> \endverbatim 149*> 150*> \param[in] LDA 151*> \verbatim 152*> LDA is INTEGER 153*> The leading dimension of the array A. LDA >= max(1,N). 154*> \endverbatim 155*> 156*> \param[in,out] AF 157*> \verbatim 158*> AF is COMPLEX*16 array, dimension (LDAF,N) 159*> If FACT = 'F', then AF is an input argument and on entry 160*> contains the triangular factor U or L from the Cholesky 161*> factorization A = U**H *U or A = L*L**H, in the same storage 162*> format as A. If EQUED .ne. 'N', then AF is the factored form 163*> of the equilibrated matrix diag(S)*A*diag(S). 164*> 165*> If FACT = 'N', then AF is an output argument and on exit 166*> returns the triangular factor U or L from the Cholesky 167*> factorization A = U**H *U or A = L*L**H of the original 168*> matrix A. 169*> 170*> If FACT = 'E', then AF is an output argument and on exit 171*> returns the triangular factor U or L from the Cholesky 172*> factorization A = U**H *U or A = L*L**H of the equilibrated 173*> matrix A (see the description of A for the form of the 174*> equilibrated matrix). 175*> \endverbatim 176*> 177*> \param[in] LDAF 178*> \verbatim 179*> LDAF is INTEGER 180*> The leading dimension of the array AF. LDAF >= max(1,N). 181*> \endverbatim 182*> 183*> \param[in,out] EQUED 184*> \verbatim 185*> EQUED is CHARACTER*1 186*> Specifies the form of equilibration that was done. 187*> = 'N': No equilibration (always true if FACT = 'N'). 188*> = 'Y': Equilibration was done, i.e., A has been replaced by 189*> diag(S) * A * diag(S). 190*> EQUED is an input argument if FACT = 'F'; otherwise, it is an 191*> output argument. 192*> \endverbatim 193*> 194*> \param[in,out] S 195*> \verbatim 196*> S is DOUBLE PRECISION array, dimension (N) 197*> The scale factors for A; not accessed if EQUED = 'N'. S is 198*> an input argument if FACT = 'F'; otherwise, S is an output 199*> argument. If FACT = 'F' and EQUED = 'Y', each element of S 200*> must be positive. 201*> \endverbatim 202*> 203*> \param[in,out] B 204*> \verbatim 205*> B is COMPLEX*16 array, dimension (LDB,NRHS) 206*> On entry, the N-by-NRHS righthand side matrix B. 207*> On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y', 208*> B is overwritten by diag(S) * B. 209*> \endverbatim 210*> 211*> \param[in] LDB 212*> \verbatim 213*> LDB is INTEGER 214*> The leading dimension of the array B. LDB >= max(1,N). 215*> \endverbatim 216*> 217*> \param[out] X 218*> \verbatim 219*> X is COMPLEX*16 array, dimension (LDX,NRHS) 220*> If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to 221*> the original system of equations. Note that if EQUED = 'Y', 222*> A and B are modified on exit, and the solution to the 223*> equilibrated system is inv(diag(S))*X. 224*> \endverbatim 225*> 226*> \param[in] LDX 227*> \verbatim 228*> LDX is INTEGER 229*> The leading dimension of the array X. LDX >= max(1,N). 230*> \endverbatim 231*> 232*> \param[out] RCOND 233*> \verbatim 234*> RCOND is DOUBLE PRECISION 235*> The estimate of the reciprocal condition number of the matrix 236*> A after equilibration (if done). If RCOND is less than the 237*> machine precision (in particular, if RCOND = 0), the matrix 238*> is singular to working precision. This condition is 239*> indicated by a return code of INFO > 0. 240*> \endverbatim 241*> 242*> \param[out] FERR 243*> \verbatim 244*> FERR is DOUBLE PRECISION array, dimension (NRHS) 245*> The estimated forward error bound for each solution vector 246*> X(j) (the j-th column of the solution matrix X). 247*> If XTRUE is the true solution corresponding to X(j), FERR(j) 248*> is an estimated upper bound for the magnitude of the largest 249*> element in (X(j) - XTRUE) divided by the magnitude of the 250*> largest element in X(j). The estimate is as reliable as 251*> the estimate for RCOND, and is almost always a slight 252*> overestimate of the true error. 253*> \endverbatim 254*> 255*> \param[out] BERR 256*> \verbatim 257*> BERR is DOUBLE PRECISION array, dimension (NRHS) 258*> The componentwise relative backward error of each solution 259*> vector X(j) (i.e., the smallest relative change in 260*> any element of A or B that makes X(j) an exact solution). 261*> \endverbatim 262*> 263*> \param[out] WORK 264*> \verbatim 265*> WORK is COMPLEX*16 array, dimension (2*N) 266*> \endverbatim 267*> 268*> \param[out] RWORK 269*> \verbatim 270*> RWORK is DOUBLE PRECISION array, dimension (N) 271*> \endverbatim 272*> 273*> \param[out] INFO 274*> \verbatim 275*> INFO is INTEGER 276*> = 0: successful exit 277*> < 0: if INFO = -i, the i-th argument had an illegal value 278*> > 0: if INFO = i, and i is 279*> <= N: the leading minor of order i of A is 280*> not positive definite, so the factorization 281*> could not be completed, and the solution has not 282*> been computed. RCOND = 0 is returned. 283*> = N+1: U is nonsingular, but RCOND is less than machine 284*> precision, meaning that the matrix is singular 285*> to working precision. Nevertheless, the 286*> solution and error bounds are computed because 287*> there are a number of situations where the 288*> computed solution can be more accurate than the 289*> value of RCOND would suggest. 290*> \endverbatim 291* 292* Authors: 293* ======== 294* 295*> \author Univ. of Tennessee 296*> \author Univ. of California Berkeley 297*> \author Univ. of Colorado Denver 298*> \author NAG Ltd. 299* 300*> \ingroup complex16POsolve 301* 302* ===================================================================== 303 SUBROUTINE ZPOSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED, 304 $ S, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, 305 $ RWORK, INFO ) 306* 307* -- LAPACK driver routine -- 308* -- LAPACK is a software package provided by Univ. of Tennessee, -- 309* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 310* 311* .. Scalar Arguments .. 312 CHARACTER EQUED, FACT, UPLO 313 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS 314 DOUBLE PRECISION RCOND 315* .. 316* .. Array Arguments .. 317 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * ), S( * ) 318 COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 319 $ WORK( * ), X( LDX, * ) 320* .. 321* 322* ===================================================================== 323* 324* .. Parameters .. 325 DOUBLE PRECISION ZERO, ONE 326 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 327* .. 328* .. Local Scalars .. 329 LOGICAL EQUIL, NOFACT, RCEQU 330 INTEGER I, INFEQU, J 331 DOUBLE PRECISION AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM 332* .. 333* .. External Functions .. 334 LOGICAL LSAME 335 DOUBLE PRECISION DLAMCH, ZLANHE 336 EXTERNAL LSAME, DLAMCH, ZLANHE 337* .. 338* .. External Subroutines .. 339 EXTERNAL XERBLA, ZLACPY, ZLAQHE, ZPOCON, ZPOEQU, ZPORFS, 340 $ ZPOTRF, ZPOTRS 341* .. 342* .. Intrinsic Functions .. 343 INTRINSIC MAX, MIN 344* .. 345* .. Executable Statements .. 346* 347 INFO = 0 348 NOFACT = LSAME( FACT, 'N' ) 349 EQUIL = LSAME( FACT, 'E' ) 350 IF( NOFACT .OR. EQUIL ) THEN 351 EQUED = 'N' 352 RCEQU = .FALSE. 353 ELSE 354 RCEQU = LSAME( EQUED, 'Y' ) 355 SMLNUM = DLAMCH( 'Safe minimum' ) 356 BIGNUM = ONE / SMLNUM 357 END IF 358* 359* Test the input parameters. 360* 361 IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) ) 362 $ THEN 363 INFO = -1 364 ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) 365 $ THEN 366 INFO = -2 367 ELSE IF( N.LT.0 ) THEN 368 INFO = -3 369 ELSE IF( NRHS.LT.0 ) THEN 370 INFO = -4 371 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 372 INFO = -6 373 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN 374 INFO = -8 375 ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT. 376 $ ( RCEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN 377 INFO = -9 378 ELSE 379 IF( RCEQU ) THEN 380 SMIN = BIGNUM 381 SMAX = ZERO 382 DO 10 J = 1, N 383 SMIN = MIN( SMIN, S( J ) ) 384 SMAX = MAX( SMAX, S( J ) ) 385 10 CONTINUE 386 IF( SMIN.LE.ZERO ) THEN 387 INFO = -10 388 ELSE IF( N.GT.0 ) THEN 389 SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM ) 390 ELSE 391 SCOND = ONE 392 END IF 393 END IF 394 IF( INFO.EQ.0 ) THEN 395 IF( LDB.LT.MAX( 1, N ) ) THEN 396 INFO = -12 397 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 398 INFO = -14 399 END IF 400 END IF 401 END IF 402* 403 IF( INFO.NE.0 ) THEN 404 CALL XERBLA( 'ZPOSVX', -INFO ) 405 RETURN 406 END IF 407* 408 IF( EQUIL ) THEN 409* 410* Compute row and column scalings to equilibrate the matrix A. 411* 412 CALL ZPOEQU( N, A, LDA, S, SCOND, AMAX, INFEQU ) 413 IF( INFEQU.EQ.0 ) THEN 414* 415* Equilibrate the matrix. 416* 417 CALL ZLAQHE( UPLO, N, A, LDA, S, SCOND, AMAX, EQUED ) 418 RCEQU = LSAME( EQUED, 'Y' ) 419 END IF 420 END IF 421* 422* Scale the right hand side. 423* 424 IF( RCEQU ) THEN 425 DO 30 J = 1, NRHS 426 DO 20 I = 1, N 427 B( I, J ) = S( I )*B( I, J ) 428 20 CONTINUE 429 30 CONTINUE 430 END IF 431* 432 IF( NOFACT .OR. EQUIL ) THEN 433* 434* Compute the Cholesky factorization A = U**H *U or A = L*L**H. 435* 436 CALL ZLACPY( UPLO, N, N, A, LDA, AF, LDAF ) 437 CALL ZPOTRF( UPLO, N, AF, LDAF, INFO ) 438* 439* Return if INFO is non-zero. 440* 441 IF( INFO.GT.0 )THEN 442 RCOND = ZERO 443 RETURN 444 END IF 445 END IF 446* 447* Compute the norm of the matrix A. 448* 449 ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK ) 450* 451* Compute the reciprocal of the condition number of A. 452* 453 CALL ZPOCON( UPLO, N, AF, LDAF, ANORM, RCOND, WORK, RWORK, INFO ) 454* 455* Compute the solution matrix X. 456* 457 CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 458 CALL ZPOTRS( UPLO, N, NRHS, AF, LDAF, X, LDX, INFO ) 459* 460* Use iterative refinement to improve the computed solution and 461* compute error bounds and backward error estimates for it. 462* 463 CALL ZPORFS( UPLO, N, NRHS, A, LDA, AF, LDAF, B, LDB, X, LDX, 464 $ FERR, BERR, WORK, RWORK, INFO ) 465* 466* Transform the solution matrix X to a solution of the original 467* system. 468* 469 IF( RCEQU ) THEN 470 DO 50 J = 1, NRHS 471 DO 40 I = 1, N 472 X( I, J ) = S( I )*X( I, J ) 473 40 CONTINUE 474 50 CONTINUE 475 DO 60 J = 1, NRHS 476 FERR( J ) = FERR( J ) / SCOND 477 60 CONTINUE 478 END IF 479* 480* Set INFO = N+1 if the matrix is singular to working precision. 481* 482 IF( RCOND.LT.DLAMCH( 'Epsilon' ) ) 483 $ INFO = N + 1 484* 485 RETURN 486* 487* End of ZPOSVX 488* 489 END 490