![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DPPSVX( FACT, UPLO, N, NRHS, AP, AFP, EQUED, S, B, LDB, 2: $ X, LDX, RCOND, FERR, BERR, WORK, IWORK, INFO ) 3: * 4: * -- LAPACK driver routine (version 3.2) -- 5: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 7: * November 2006 8: * 9: * .. Scalar Arguments .. 10: CHARACTER EQUED, FACT, UPLO 11: INTEGER INFO, LDB, LDX, N, NRHS 12: DOUBLE PRECISION RCOND 13: * .. 14: * .. Array Arguments .. 15: INTEGER IWORK( * ) 16: DOUBLE PRECISION AFP( * ), AP( * ), B( LDB, * ), BERR( * ), 17: $ FERR( * ), S( * ), WORK( * ), X( LDX, * ) 18: * .. 19: * 20: * Purpose 21: * ======= 22: * 23: * DPPSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to 24: * compute the solution to a real system of linear equations 25: * A * X = B, 26: * where A is an N-by-N symmetric positive definite matrix stored in 27: * packed format and X and B are N-by-NRHS matrices. 28: * 29: * Error bounds on the solution and a condition estimate are also 30: * provided. 31: * 32: * Description 33: * =========== 34: * 35: * The following steps are performed: 36: * 37: * 1. If FACT = 'E', real scaling factors are computed to equilibrate 38: * the system: 39: * diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B 40: * Whether or not the system will be equilibrated depends on the 41: * scaling of the matrix A, but if equilibration is used, A is 42: * overwritten by diag(S)*A*diag(S) and B by diag(S)*B. 43: * 44: * 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to 45: * factor the matrix A (after equilibration if FACT = 'E') as 46: * A = U**T* U, if UPLO = 'U', or 47: * A = L * L**T, if UPLO = 'L', 48: * where U is an upper triangular matrix and L is a lower triangular 49: * matrix. 50: * 51: * 3. If the leading i-by-i principal minor is not positive definite, 52: * then the routine returns with INFO = i. Otherwise, the factored 53: * form of A is used to estimate the condition number of the matrix 54: * A. If the reciprocal of the condition number is less than machine 55: * precision, INFO = N+1 is returned as a warning, but the routine 56: * still goes on to solve for X and compute error bounds as 57: * described below. 58: * 59: * 4. The system of equations is solved for X using the factored form 60: * of A. 61: * 62: * 5. Iterative refinement is applied to improve the computed solution 63: * matrix and calculate error bounds and backward error estimates 64: * for it. 65: * 66: * 6. If equilibration was used, the matrix X is premultiplied by 67: * diag(S) so that it solves the original system before 68: * equilibration. 69: * 70: * Arguments 71: * ========= 72: * 73: * FACT (input) CHARACTER*1 74: * Specifies whether or not the factored form of the matrix A is 75: * supplied on entry, and if not, whether the matrix A should be 76: * equilibrated before it is factored. 77: * = 'F': On entry, AFP contains the factored form of A. 78: * If EQUED = 'Y', the matrix A has been equilibrated 79: * with scaling factors given by S. AP and AFP will not 80: * be modified. 81: * = 'N': The matrix A will be copied to AFP and factored. 82: * = 'E': The matrix A will be equilibrated if necessary, then 83: * copied to AFP and factored. 84: * 85: * UPLO (input) CHARACTER*1 86: * = 'U': Upper triangle of A is stored; 87: * = 'L': Lower triangle of A is stored. 88: * 89: * N (input) INTEGER 90: * The number of linear equations, i.e., the order of the 91: * matrix A. N >= 0. 92: * 93: * NRHS (input) INTEGER 94: * The number of right hand sides, i.e., the number of columns 95: * of the matrices B and X. NRHS >= 0. 96: * 97: * AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2) 98: * On entry, the upper or lower triangle of the symmetric matrix 99: * A, packed columnwise in a linear array, except if FACT = 'F' 100: * and EQUED = 'Y', then A must contain the equilibrated matrix 101: * diag(S)*A*diag(S). The j-th column of A is stored in the 102: * array AP as follows: 103: * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 104: * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. 105: * See below for further details. A is not modified if 106: * FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit. 107: * 108: * On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by 109: * diag(S)*A*diag(S). 110: * 111: * AFP (input or output) DOUBLE PRECISION array, dimension 112: * (N*(N+1)/2) 113: * If FACT = 'F', then AFP is an input argument and on entry 114: * contains the triangular factor U or L from the Cholesky 115: * factorization A = U'*U or A = L*L', in the same storage 116: * format as A. If EQUED .ne. 'N', then AFP is the factored 117: * form of the equilibrated matrix A. 118: * 119: * If FACT = 'N', then AFP is an output argument and on exit 120: * returns the triangular factor U or L from the Cholesky 121: * factorization A = U'*U or A = L*L' of the original matrix A. 122: * 123: * If FACT = 'E', then AFP is an output argument and on exit 124: * returns the triangular factor U or L from the Cholesky 125: * factorization A = U'*U or A = L*L' of the equilibrated 126: * matrix A (see the description of AP for the form of the 127: * equilibrated matrix). 128: * 129: * EQUED (input or output) CHARACTER*1 130: * Specifies the form of equilibration that was done. 131: * = 'N': No equilibration (always true if FACT = 'N'). 132: * = 'Y': Equilibration was done, i.e., A has been replaced by 133: * diag(S) * A * diag(S). 134: * EQUED is an input argument if FACT = 'F'; otherwise, it is an 135: * output argument. 136: * 137: * S (input or output) DOUBLE PRECISION array, dimension (N) 138: * The scale factors for A; not accessed if EQUED = 'N'. S is 139: * an input argument if FACT = 'F'; otherwise, S is an output 140: * argument. If FACT = 'F' and EQUED = 'Y', each element of S 141: * must be positive. 142: * 143: * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) 144: * On entry, the N-by-NRHS right hand side matrix B. 145: * On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y', 146: * B is overwritten by diag(S) * B. 147: * 148: * LDB (input) INTEGER 149: * The leading dimension of the array B. LDB >= max(1,N). 150: * 151: * X (output) DOUBLE PRECISION array, dimension (LDX,NRHS) 152: * If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to 153: * the original system of equations. Note that if EQUED = 'Y', 154: * A and B are modified on exit, and the solution to the 155: * equilibrated system is inv(diag(S))*X. 156: * 157: * LDX (input) INTEGER 158: * The leading dimension of the array X. LDX >= max(1,N). 159: * 160: * RCOND (output) DOUBLE PRECISION 161: * The estimate of the reciprocal condition number of the matrix 162: * A after equilibration (if done). If RCOND is less than the 163: * machine precision (in particular, if RCOND = 0), the matrix 164: * is singular to working precision. This condition is 165: * indicated by a return code of INFO > 0. 166: * 167: * FERR (output) DOUBLE PRECISION array, dimension (NRHS) 168: * The estimated forward error bound for each solution vector 169: * X(j) (the j-th column of the solution matrix X). 170: * If XTRUE is the true solution corresponding to X(j), FERR(j) 171: * is an estimated upper bound for the magnitude of the largest 172: * element in (X(j) - XTRUE) divided by the magnitude of the 173: * largest element in X(j). The estimate is as reliable as 174: * the estimate for RCOND, and is almost always a slight 175: * overestimate of the true error. 176: * 177: * BERR (output) DOUBLE PRECISION array, dimension (NRHS) 178: * The componentwise relative backward error of each solution 179: * vector X(j) (i.e., the smallest relative change in 180: * any element of A or B that makes X(j) an exact solution). 181: * 182: * WORK (workspace) DOUBLE PRECISION array, dimension (3*N) 183: * 184: * IWORK (workspace) INTEGER array, dimension (N) 185: * 186: * INFO (output) INTEGER 187: * = 0: successful exit 188: * < 0: if INFO = -i, the i-th argument had an illegal value 189: * > 0: if INFO = i, and i is 190: * <= N: the leading minor of order i of A is 191: * not positive definite, so the factorization 192: * could not be completed, and the solution has not 193: * been computed. RCOND = 0 is returned. 194: * = N+1: U is nonsingular, but RCOND is less than machine 195: * precision, meaning that the matrix is singular 196: * to working precision. Nevertheless, the 197: * solution and error bounds are computed because 198: * there are a number of situations where the 199: * computed solution can be more accurate than the 200: * value of RCOND would suggest. 201: * 202: * Further Details 203: * =============== 204: * 205: * The packed storage scheme is illustrated by the following example 206: * when N = 4, UPLO = 'U': 207: * 208: * Two-dimensional storage of the symmetric matrix A: 209: * 210: * a11 a12 a13 a14 211: * a22 a23 a24 212: * a33 a34 (aij = conjg(aji)) 213: * a44 214: * 215: * Packed storage of the upper triangle of A: 216: * 217: * AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ] 218: * 219: * ===================================================================== 220: * 221: * .. Parameters .. 222: DOUBLE PRECISION ZERO, ONE 223: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 224: * .. 225: * .. Local Scalars .. 226: LOGICAL EQUIL, NOFACT, RCEQU 227: INTEGER I, INFEQU, J 228: DOUBLE PRECISION AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM 229: * .. 230: * .. External Functions .. 231: LOGICAL LSAME 232: DOUBLE PRECISION DLAMCH, DLANSP 233: EXTERNAL LSAME, DLAMCH, DLANSP 234: * .. 235: * .. External Subroutines .. 236: EXTERNAL DCOPY, DLACPY, DLAQSP, DPPCON, DPPEQU, DPPRFS, 237: $ DPPTRF, DPPTRS, XERBLA 238: * .. 239: * .. Intrinsic Functions .. 240: INTRINSIC MAX, MIN 241: * .. 242: * .. Executable Statements .. 243: * 244: INFO = 0 245: NOFACT = LSAME( FACT, 'N' ) 246: EQUIL = LSAME( FACT, 'E' ) 247: IF( NOFACT .OR. EQUIL ) THEN 248: EQUED = 'N' 249: RCEQU = .FALSE. 250: ELSE 251: RCEQU = LSAME( EQUED, 'Y' ) 252: SMLNUM = DLAMCH( 'Safe minimum' ) 253: BIGNUM = ONE / SMLNUM 254: END IF 255: * 256: * Test the input parameters. 257: * 258: IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) ) 259: $ THEN 260: INFO = -1 261: ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) 262: $ THEN 263: INFO = -2 264: ELSE IF( N.LT.0 ) THEN 265: INFO = -3 266: ELSE IF( NRHS.LT.0 ) THEN 267: INFO = -4 268: ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT. 269: $ ( RCEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN 270: INFO = -7 271: ELSE 272: IF( RCEQU ) THEN 273: SMIN = BIGNUM 274: SMAX = ZERO 275: DO 10 J = 1, N 276: SMIN = MIN( SMIN, S( J ) ) 277: SMAX = MAX( SMAX, S( J ) ) 278: 10 CONTINUE 279: IF( SMIN.LE.ZERO ) THEN 280: INFO = -8 281: ELSE IF( N.GT.0 ) THEN 282: SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM ) 283: ELSE 284: SCOND = ONE 285: END IF 286: END IF 287: IF( INFO.EQ.0 ) THEN 288: IF( LDB.LT.MAX( 1, N ) ) THEN 289: INFO = -10 290: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 291: INFO = -12 292: END IF 293: END IF 294: END IF 295: * 296: IF( INFO.NE.0 ) THEN 297: CALL XERBLA( 'DPPSVX', -INFO ) 298: RETURN 299: END IF 300: * 301: IF( EQUIL ) THEN 302: * 303: * Compute row and column scalings to equilibrate the matrix A. 304: * 305: CALL DPPEQU( UPLO, N, AP, S, SCOND, AMAX, INFEQU ) 306: IF( INFEQU.EQ.0 ) THEN 307: * 308: * Equilibrate the matrix. 309: * 310: CALL DLAQSP( UPLO, N, AP, S, SCOND, AMAX, EQUED ) 311: RCEQU = LSAME( EQUED, 'Y' ) 312: END IF 313: END IF 314: * 315: * Scale the right-hand side. 316: * 317: IF( RCEQU ) THEN 318: DO 30 J = 1, NRHS 319: DO 20 I = 1, N 320: B( I, J ) = S( I )*B( I, J ) 321: 20 CONTINUE 322: 30 CONTINUE 323: END IF 324: * 325: IF( NOFACT .OR. EQUIL ) THEN 326: * 327: * Compute the Cholesky factorization A = U'*U or A = L*L'. 328: * 329: CALL DCOPY( N*( N+1 ) / 2, AP, 1, AFP, 1 ) 330: CALL DPPTRF( UPLO, N, AFP, INFO ) 331: * 332: * Return if INFO is non-zero. 333: * 334: IF( INFO.GT.0 )THEN 335: RCOND = ZERO 336: RETURN 337: END IF 338: END IF 339: * 340: * Compute the norm of the matrix A. 341: * 342: ANORM = DLANSP( 'I', UPLO, N, AP, WORK ) 343: * 344: * Compute the reciprocal of the condition number of A. 345: * 346: CALL DPPCON( UPLO, N, AFP, ANORM, RCOND, WORK, IWORK, INFO ) 347: * 348: * Compute the solution matrix X. 349: * 350: CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 351: CALL DPPTRS( UPLO, N, NRHS, AFP, X, LDX, INFO ) 352: * 353: * Use iterative refinement to improve the computed solution and 354: * compute error bounds and backward error estimates for it. 355: * 356: CALL DPPRFS( UPLO, N, NRHS, AP, AFP, B, LDB, X, LDX, FERR, BERR, 357: $ WORK, IWORK, INFO ) 358: * 359: * Transform the solution matrix X to a solution of the original 360: * system. 361: * 362: IF( RCEQU ) THEN 363: DO 50 J = 1, NRHS 364: DO 40 I = 1, N 365: X( I, J ) = S( I )*X( I, J ) 366: 40 CONTINUE 367: 50 CONTINUE 368: DO 60 J = 1, NRHS 369: FERR( J ) = FERR( J ) / SCOND 370: 60 CONTINUE 371: END IF 372: * 373: * Set INFO = N+1 if the matrix is singular to working precision. 374: * 375: IF( RCOND.LT.DLAMCH( 'Epsilon' ) ) 376: $ INFO = N + 1 377: * 378: RETURN 379: * 380: * End of DPPSVX 381: * 382: END