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