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