![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DSYSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, 2: $ LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK, 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 FACT, UPLO 12: INTEGER INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS 13: DOUBLE PRECISION RCOND 14: * .. 15: * .. Array Arguments .. 16: INTEGER IPIV( * ), IWORK( * ) 17: DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 18: $ BERR( * ), FERR( * ), WORK( * ), X( LDX, * ) 19: * .. 20: * 21: * Purpose 22: * ======= 23: * 24: * DSYSVX uses the diagonal pivoting factorization to compute the 25: * solution to a real system of linear equations A * X = B, 26: * where A is an N-by-N symmetric matrix and X and B are N-by-NRHS 27: * 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 = 'N', the diagonal pivoting method is used to factor A. 38: * The form of the factorization is 39: * A = U * D * U**T, if UPLO = 'U', or 40: * A = L * D * L**T, if UPLO = 'L', 41: * where U (or L) is a product of permutation and unit upper (lower) 42: * triangular matrices, and D is symmetric and block diagonal with 43: * 1-by-1 and 2-by-2 diagonal blocks. 44: * 45: * 2. If some D(i,i)=0, so that D is exactly singular, then the routine 46: * returns with INFO = i. Otherwise, the factored form of A is used 47: * to estimate the condition number of the matrix A. If the 48: * reciprocal of the condition number is less than machine precision, 49: * INFO = N+1 is returned as a warning, but the routine still goes on 50: * to solve for X and compute error bounds as described below. 51: * 52: * 3. The system of equations is solved for X using the factored form 53: * of A. 54: * 55: * 4. Iterative refinement is applied to improve the computed solution 56: * matrix and calculate error bounds and backward error estimates 57: * for it. 58: * 59: * Arguments 60: * ========= 61: * 62: * FACT (input) CHARACTER*1 63: * Specifies whether or not the factored form of A has been 64: * supplied on entry. 65: * = 'F': On entry, AF and IPIV contain the factored form of 66: * A. AF and IPIV will not be modified. 67: * = 'N': The matrix A will be copied to AF and factored. 68: * 69: * UPLO (input) CHARACTER*1 70: * = 'U': Upper triangle of A is stored; 71: * = 'L': Lower triangle of A is stored. 72: * 73: * N (input) INTEGER 74: * The number of linear equations, i.e., the order of the 75: * matrix A. N >= 0. 76: * 77: * NRHS (input) INTEGER 78: * The number of right hand sides, i.e., the number of columns 79: * of the matrices B and X. NRHS >= 0. 80: * 81: * A (input) DOUBLE PRECISION array, dimension (LDA,N) 82: * The symmetric matrix A. If UPLO = 'U', the leading N-by-N 83: * upper triangular part of A contains the upper triangular part 84: * of the matrix A, and the strictly lower triangular part of A 85: * is not referenced. If UPLO = 'L', the leading N-by-N lower 86: * triangular part of A contains the lower triangular part of 87: * the matrix A, and the strictly upper triangular part of A is 88: * not referenced. 89: * 90: * LDA (input) INTEGER 91: * The leading dimension of the array A. LDA >= max(1,N). 92: * 93: * AF (input or output) DOUBLE PRECISION array, dimension (LDAF,N) 94: * If FACT = 'F', then AF is an input argument and on entry 95: * contains the block diagonal matrix D and the multipliers used 96: * to obtain the factor U or L from the factorization 97: * A = U*D*U**T or A = L*D*L**T as computed by DSYTRF. 98: * 99: * If FACT = 'N', then AF is an output argument and on exit 100: * returns the block diagonal matrix D and the multipliers used 101: * to obtain the factor U or L from the factorization 102: * A = U*D*U**T or A = L*D*L**T. 103: * 104: * LDAF (input) INTEGER 105: * The leading dimension of the array AF. LDAF >= max(1,N). 106: * 107: * IPIV (input or output) INTEGER array, dimension (N) 108: * If FACT = 'F', then IPIV is an input argument and on entry 109: * contains details of the interchanges and the block structure 110: * of D, as determined by DSYTRF. 111: * If IPIV(k) > 0, then rows and columns k and IPIV(k) were 112: * interchanged and D(k,k) is a 1-by-1 diagonal block. 113: * If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and 114: * columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) 115: * is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = 116: * IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were 117: * interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. 118: * 119: * If FACT = 'N', then IPIV is an output argument and on exit 120: * contains details of the interchanges and the block structure 121: * of D, as determined by DSYTRF. 122: * 123: * B (input) DOUBLE PRECISION array, dimension (LDB,NRHS) 124: * The N-by-NRHS right hand side matrix B. 125: * 126: * LDB (input) INTEGER 127: * The leading dimension of the array B. LDB >= max(1,N). 128: * 129: * X (output) DOUBLE PRECISION array, dimension (LDX,NRHS) 130: * If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X. 131: * 132: * LDX (input) INTEGER 133: * The leading dimension of the array X. LDX >= max(1,N). 134: * 135: * RCOND (output) DOUBLE PRECISION 136: * The estimate of the reciprocal condition number of the matrix 137: * A. If RCOND is less than the machine precision (in 138: * particular, if RCOND = 0), the matrix is singular to working 139: * precision. This condition is indicated by a return code of 140: * INFO > 0. 141: * 142: * FERR (output) DOUBLE PRECISION array, dimension (NRHS) 143: * The estimated forward error bound for each solution vector 144: * X(j) (the j-th column of the solution matrix X). 145: * If XTRUE is the true solution corresponding to X(j), FERR(j) 146: * is an estimated upper bound for the magnitude of the largest 147: * element in (X(j) - XTRUE) divided by the magnitude of the 148: * largest element in X(j). The estimate is as reliable as 149: * the estimate for RCOND, and is almost always a slight 150: * overestimate of the true error. 151: * 152: * BERR (output) DOUBLE PRECISION array, dimension (NRHS) 153: * The componentwise relative backward error of each solution 154: * vector X(j) (i.e., the smallest relative change in 155: * any element of A or B that makes X(j) an exact solution). 156: * 157: * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 158: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 159: * 160: * LWORK (input) INTEGER 161: * The length of WORK. LWORK >= max(1,3*N), and for best 162: * performance, when FACT = 'N', LWORK >= max(1,3*N,N*NB), where 163: * NB is the optimal blocksize for DSYTRF. 164: * 165: * If LWORK = -1, then a workspace query is assumed; the routine 166: * only calculates the optimal size of the WORK array, returns 167: * this value as the first entry of the WORK array, and no error 168: * message related to LWORK is issued by XERBLA. 169: * 170: * IWORK (workspace) INTEGER array, dimension (N) 171: * 172: * INFO (output) INTEGER 173: * = 0: successful exit 174: * < 0: if INFO = -i, the i-th argument had an illegal value 175: * > 0: if INFO = i, and i is 176: * <= N: D(i,i) is exactly zero. The factorization 177: * has been completed but the factor D is exactly 178: * singular, so the solution and error bounds could 179: * not be computed. RCOND = 0 is returned. 180: * = N+1: D is nonsingular, but RCOND is less than machine 181: * precision, meaning that the matrix is singular 182: * to working precision. Nevertheless, the 183: * solution and error bounds are computed because 184: * there are a number of situations where the 185: * computed solution can be more accurate than the 186: * value of RCOND would suggest. 187: * 188: * ===================================================================== 189: * 190: * .. Parameters .. 191: DOUBLE PRECISION ZERO 192: PARAMETER ( ZERO = 0.0D+0 ) 193: * .. 194: * .. Local Scalars .. 195: LOGICAL LQUERY, NOFACT 196: INTEGER LWKOPT, NB 197: DOUBLE PRECISION ANORM 198: * .. 199: * .. External Functions .. 200: LOGICAL LSAME 201: INTEGER ILAENV 202: DOUBLE PRECISION DLAMCH, DLANSY 203: EXTERNAL LSAME, ILAENV, DLAMCH, DLANSY 204: * .. 205: * .. External Subroutines .. 206: EXTERNAL DLACPY, DSYCON, DSYRFS, DSYTRF, DSYTRS, XERBLA 207: * .. 208: * .. Intrinsic Functions .. 209: INTRINSIC MAX 210: * .. 211: * .. Executable Statements .. 212: * 213: * Test the input parameters. 214: * 215: INFO = 0 216: NOFACT = LSAME( FACT, 'N' ) 217: LQUERY = ( LWORK.EQ.-1 ) 218: IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN 219: INFO = -1 220: ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) 221: $ THEN 222: INFO = -2 223: ELSE IF( N.LT.0 ) THEN 224: INFO = -3 225: ELSE IF( NRHS.LT.0 ) THEN 226: INFO = -4 227: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 228: INFO = -6 229: ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN 230: INFO = -8 231: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 232: INFO = -11 233: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 234: INFO = -13 235: ELSE IF( LWORK.LT.MAX( 1, 3*N ) .AND. .NOT.LQUERY ) THEN 236: INFO = -18 237: END IF 238: * 239: IF( INFO.EQ.0 ) THEN 240: LWKOPT = MAX( 1, 3*N ) 241: IF( NOFACT ) THEN 242: NB = ILAENV( 1, 'DSYTRF', UPLO, N, -1, -1, -1 ) 243: LWKOPT = MAX( LWKOPT, N*NB ) 244: END IF 245: WORK( 1 ) = LWKOPT 246: END IF 247: * 248: IF( INFO.NE.0 ) THEN 249: CALL XERBLA( 'DSYSVX', -INFO ) 250: RETURN 251: ELSE IF( LQUERY ) THEN 252: RETURN 253: END IF 254: * 255: IF( NOFACT ) THEN 256: * 257: * Compute the factorization A = U*D*U' or A = L*D*L'. 258: * 259: CALL DLACPY( UPLO, N, N, A, LDA, AF, LDAF ) 260: CALL DSYTRF( UPLO, N, AF, LDAF, IPIV, WORK, LWORK, INFO ) 261: * 262: * Return if INFO is non-zero. 263: * 264: IF( INFO.GT.0 )THEN 265: RCOND = ZERO 266: RETURN 267: END IF 268: END IF 269: * 270: * Compute the norm of the matrix A. 271: * 272: ANORM = DLANSY( 'I', UPLO, N, A, LDA, WORK ) 273: * 274: * Compute the reciprocal of the condition number of A. 275: * 276: CALL DSYCON( UPLO, N, AF, LDAF, IPIV, ANORM, RCOND, WORK, IWORK, 277: $ INFO ) 278: * 279: * Compute the solution vectors X. 280: * 281: CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 282: CALL DSYTRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO ) 283: * 284: * Use iterative refinement to improve the computed solutions and 285: * compute error bounds and backward error estimates for them. 286: * 287: CALL DSYRFS( UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, 288: $ LDX, FERR, BERR, WORK, IWORK, INFO ) 289: * 290: * Set INFO = N+1 if the matrix is singular to working precision. 291: * 292: IF( RCOND.LT.DLAMCH( 'Epsilon' ) ) 293: $ INFO = N + 1 294: * 295: WORK( 1 ) = LWKOPT 296: * 297: RETURN 298: * 299: * End of DSYSVX 300: * 301: END