1*> \brief <b> CPOSVX 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 CPOSVX + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cposvx.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cposvx.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cposvx.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE CPOSVX( 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* REAL RCOND 29* .. 30* .. Array Arguments .. 31* REAL BERR( * ), FERR( * ), RWORK( * ), S( * ) 32* COMPLEX A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 33* $ WORK( * ), X( LDX, * ) 34* .. 35* 36* 37*> \par Purpose: 38* ============= 39*> 40*> \verbatim 41*> 42*> CPOSVX 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 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 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 REAL 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 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 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 REAL 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 REAL 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 REAL 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 array, dimension (2*N) 266*> \endverbatim 267*> 268*> \param[out] RWORK 269*> \verbatim 270*> RWORK is REAL 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*> \date April 2012 301* 302*> \ingroup complexPOsolve 303* 304* ===================================================================== 305 SUBROUTINE CPOSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED, 306 $ S, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, 307 $ RWORK, INFO ) 308* 309* -- LAPACK driver routine (version 3.4.1) -- 310* -- LAPACK is a software package provided by Univ. of Tennessee, -- 311* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 312* April 2012 313* 314* .. Scalar Arguments .. 315 CHARACTER EQUED, FACT, UPLO 316 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS 317 REAL RCOND 318* .. 319* .. Array Arguments .. 320 REAL BERR( * ), FERR( * ), RWORK( * ), S( * ) 321 COMPLEX A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 322 $ WORK( * ), X( LDX, * ) 323* .. 324* 325* ===================================================================== 326* 327* .. Parameters .. 328 REAL ZERO, ONE 329 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 330* .. 331* .. Local Scalars .. 332 LOGICAL EQUIL, NOFACT, RCEQU 333 INTEGER I, INFEQU, J 334 REAL AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM 335* .. 336* .. External Functions .. 337 LOGICAL LSAME 338 REAL CLANHE, SLAMCH 339 EXTERNAL LSAME, CLANHE, SLAMCH 340* .. 341* .. External Subroutines .. 342 EXTERNAL CLACPY, CLAQHE, CPOCON, CPOEQU, CPORFS, CPOTRF, 343 $ CPOTRS, XERBLA 344* .. 345* .. Intrinsic Functions .. 346 INTRINSIC MAX, MIN 347* .. 348* .. Executable Statements .. 349* 350 INFO = 0 351 NOFACT = LSAME( FACT, 'N' ) 352 EQUIL = LSAME( FACT, 'E' ) 353 IF( NOFACT .OR. EQUIL ) THEN 354 EQUED = 'N' 355 RCEQU = .FALSE. 356 ELSE 357 RCEQU = LSAME( EQUED, 'Y' ) 358 SMLNUM = SLAMCH( 'Safe minimum' ) 359 BIGNUM = ONE / SMLNUM 360 END IF 361* 362* Test the input parameters. 363* 364 IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) ) 365 $ THEN 366 INFO = -1 367 ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) 368 $ THEN 369 INFO = -2 370 ELSE IF( N.LT.0 ) THEN 371 INFO = -3 372 ELSE IF( NRHS.LT.0 ) THEN 373 INFO = -4 374 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 375 INFO = -6 376 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN 377 INFO = -8 378 ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT. 379 $ ( RCEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN 380 INFO = -9 381 ELSE 382 IF( RCEQU ) THEN 383 SMIN = BIGNUM 384 SMAX = ZERO 385 DO 10 J = 1, N 386 SMIN = MIN( SMIN, S( J ) ) 387 SMAX = MAX( SMAX, S( J ) ) 388 10 CONTINUE 389 IF( SMIN.LE.ZERO ) THEN 390 INFO = -10 391 ELSE IF( N.GT.0 ) THEN 392 SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM ) 393 ELSE 394 SCOND = ONE 395 END IF 396 END IF 397 IF( INFO.EQ.0 ) THEN 398 IF( LDB.LT.MAX( 1, N ) ) THEN 399 INFO = -12 400 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 401 INFO = -14 402 END IF 403 END IF 404 END IF 405* 406 IF( INFO.NE.0 ) THEN 407 CALL XERBLA( 'CPOSVX', -INFO ) 408 RETURN 409 END IF 410* 411 IF( EQUIL ) THEN 412* 413* Compute row and column scalings to equilibrate the matrix A. 414* 415 CALL CPOEQU( N, A, LDA, S, SCOND, AMAX, INFEQU ) 416 IF( INFEQU.EQ.0 ) THEN 417* 418* Equilibrate the matrix. 419* 420 CALL CLAQHE( UPLO, N, A, LDA, S, SCOND, AMAX, EQUED ) 421 RCEQU = LSAME( EQUED, 'Y' ) 422 END IF 423 END IF 424* 425* Scale the right hand side. 426* 427 IF( RCEQU ) THEN 428 DO 30 J = 1, NRHS 429 DO 20 I = 1, N 430 B( I, J ) = S( I )*B( I, J ) 431 20 CONTINUE 432 30 CONTINUE 433 END IF 434* 435 IF( NOFACT .OR. EQUIL ) THEN 436* 437* Compute the Cholesky factorization A = U**H *U or A = L*L**H. 438* 439 CALL CLACPY( UPLO, N, N, A, LDA, AF, LDAF ) 440 CALL CPOTRF( UPLO, N, AF, LDAF, INFO ) 441* 442* Return if INFO is non-zero. 443* 444 IF( INFO.GT.0 )THEN 445 RCOND = ZERO 446 RETURN 447 END IF 448 END IF 449* 450* Compute the norm of the matrix A. 451* 452 ANORM = CLANHE( '1', UPLO, N, A, LDA, RWORK ) 453* 454* Compute the reciprocal of the condition number of A. 455* 456 CALL CPOCON( UPLO, N, AF, LDAF, ANORM, RCOND, WORK, RWORK, INFO ) 457* 458* Compute the solution matrix X. 459* 460 CALL CLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 461 CALL CPOTRS( UPLO, N, NRHS, AF, LDAF, X, LDX, INFO ) 462* 463* Use iterative refinement to improve the computed solution and 464* compute error bounds and backward error estimates for it. 465* 466 CALL CPORFS( UPLO, N, NRHS, A, LDA, AF, LDAF, B, LDB, X, LDX, 467 $ FERR, BERR, WORK, RWORK, INFO ) 468* 469* Transform the solution matrix X to a solution of the original 470* system. 471* 472 IF( RCEQU ) THEN 473 DO 50 J = 1, NRHS 474 DO 40 I = 1, N 475 X( I, J ) = S( I )*X( I, J ) 476 40 CONTINUE 477 50 CONTINUE 478 DO 60 J = 1, NRHS 479 FERR( J ) = FERR( J ) / SCOND 480 60 CONTINUE 481 END IF 482* 483* Set INFO = N+1 if the matrix is singular to working precision. 484* 485 IF( RCOND.LT.SLAMCH( 'Epsilon' ) ) 486 $ INFO = N + 1 487* 488 RETURN 489* 490* End of CPOSVX 491* 492 END 493