![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DSPOSV( UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, WORK, 2: + SWORK, ITER, INFO ) 3: * 4: * -- LAPACK PROTOTYPE driver routine (version 3.3.0) -- 5: * Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd.. 6: * November 2010 7: * 8: * .. 9: * .. Scalar Arguments .. 10: CHARACTER UPLO 11: INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS 12: * .. 13: * .. Array Arguments .. 14: REAL SWORK( * ) 15: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( N, * ), 16: + X( LDX, * ) 17: * .. 18: * 19: * Purpose 20: * ======= 21: * 22: * DSPOSV computes the solution to a real system of linear equations 23: * A * X = B, 24: * where A is an N-by-N symmetric positive definite matrix and X and B 25: * are N-by-NRHS matrices. 26: * 27: * DSPOSV first attempts to factorize the matrix in SINGLE PRECISION 28: * and use this factorization within an iterative refinement procedure 29: * to produce a solution with DOUBLE PRECISION normwise backward error 30: * quality (see below). If the approach fails the method switches to a 31: * DOUBLE PRECISION factorization and solve. 32: * 33: * The iterative refinement is not going to be a winning strategy if 34: * the ratio SINGLE PRECISION performance over DOUBLE PRECISION 35: * performance is too small. A reasonable strategy should take the 36: * number of right-hand sides and the size of the matrix into account. 37: * This might be done with a call to ILAENV in the future. Up to now, we 38: * always try iterative refinement. 39: * 40: * The iterative refinement process is stopped if 41: * ITER > ITERMAX 42: * or for all the RHS we have: 43: * RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX 44: * where 45: * o ITER is the number of the current iteration in the iterative 46: * refinement process 47: * o RNRM is the infinity-norm of the residual 48: * o XNRM is the infinity-norm of the solution 49: * o ANRM is the infinity-operator-norm of the matrix A 50: * o EPS is the machine epsilon returned by DLAMCH('Epsilon') 51: * The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00 52: * respectively. 53: * 54: * Arguments 55: * ========= 56: * 57: * UPLO (input) CHARACTER*1 58: * = 'U': Upper triangle of A is stored; 59: * = 'L': Lower triangle of A is stored. 60: * 61: * N (input) INTEGER 62: * The number of linear equations, i.e., the order of the 63: * matrix A. N >= 0. 64: * 65: * NRHS (input) INTEGER 66: * The number of right hand sides, i.e., the number of columns 67: * of the matrix B. NRHS >= 0. 68: * 69: * A (input/output) DOUBLE PRECISION array, 70: * dimension (LDA,N) 71: * On entry, the symmetric matrix A. If UPLO = 'U', the leading 72: * N-by-N upper triangular part of A contains the upper 73: * triangular part of the matrix A, and the strictly lower 74: * triangular part of A is not referenced. If UPLO = 'L', the 75: * leading N-by-N lower triangular part of A contains the lower 76: * triangular part of the matrix A, and the strictly upper 77: * triangular part of A is not referenced. 78: * On exit, if iterative refinement has been successfully used 79: * (INFO.EQ.0 and ITER.GE.0, see description below), then A is 80: * unchanged, if double precision factorization has been used 81: * (INFO.EQ.0 and ITER.LT.0, see description below), then the 82: * array A contains the factor U or L from the Cholesky 83: * factorization A = U**T*U or A = L*L**T. 84: * 85: * 86: * LDA (input) INTEGER 87: * The leading dimension of the array A. LDA >= max(1,N). 88: * 89: * B (input) DOUBLE PRECISION array, dimension (LDB,NRHS) 90: * The N-by-NRHS right hand side matrix B. 91: * 92: * LDB (input) INTEGER 93: * The leading dimension of the array B. LDB >= max(1,N). 94: * 95: * X (output) DOUBLE PRECISION array, dimension (LDX,NRHS) 96: * If INFO = 0, the N-by-NRHS solution matrix X. 97: * 98: * LDX (input) INTEGER 99: * The leading dimension of the array X. LDX >= max(1,N). 100: * 101: * WORK (workspace) DOUBLE PRECISION array, dimension (N,NRHS) 102: * This array is used to hold the residual vectors. 103: * 104: * SWORK (workspace) REAL array, dimension (N*(N+NRHS)) 105: * This array is used to use the single precision matrix and the 106: * right-hand sides or solutions in single precision. 107: * 108: * ITER (output) INTEGER 109: * < 0: iterative refinement has failed, double precision 110: * factorization has been performed 111: * -1 : the routine fell back to full precision for 112: * implementation- or machine-specific reasons 113: * -2 : narrowing the precision induced an overflow, 114: * the routine fell back to full precision 115: * -3 : failure of SPOTRF 116: * -31: stop the iterative refinement after the 30th 117: * iterations 118: * > 0: iterative refinement has been sucessfully used. 119: * Returns the number of iterations 120: * 121: * INFO (output) INTEGER 122: * = 0: successful exit 123: * < 0: if INFO = -i, the i-th argument had an illegal value 124: * > 0: if INFO = i, the leading minor of order i of (DOUBLE 125: * PRECISION) A is not positive definite, so the 126: * factorization could not be completed, and the solution 127: * has not been computed. 128: * 129: * ========= 130: * 131: * .. Parameters .. 132: LOGICAL DOITREF 133: PARAMETER ( DOITREF = .TRUE. ) 134: * 135: INTEGER ITERMAX 136: PARAMETER ( ITERMAX = 30 ) 137: * 138: DOUBLE PRECISION BWDMAX 139: PARAMETER ( BWDMAX = 1.0E+00 ) 140: * 141: DOUBLE PRECISION NEGONE, ONE 142: PARAMETER ( NEGONE = -1.0D+0, ONE = 1.0D+0 ) 143: * 144: * .. Local Scalars .. 145: INTEGER I, IITER, PTSA, PTSX 146: DOUBLE PRECISION ANRM, CTE, EPS, RNRM, XNRM 147: * 148: * .. External Subroutines .. 149: EXTERNAL DAXPY, DSYMM, DLACPY, DLAT2S, DLAG2S, SLAG2D, 150: + SPOTRF, SPOTRS, XERBLA 151: * .. 152: * .. External Functions .. 153: INTEGER IDAMAX 154: DOUBLE PRECISION DLAMCH, DLANSY 155: LOGICAL LSAME 156: EXTERNAL IDAMAX, DLAMCH, DLANSY, LSAME 157: * .. 158: * .. Intrinsic Functions .. 159: INTRINSIC ABS, DBLE, MAX, SQRT 160: * .. 161: * .. Executable Statements .. 162: * 163: INFO = 0 164: ITER = 0 165: * 166: * Test the input parameters. 167: * 168: IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 169: INFO = -1 170: ELSE IF( N.LT.0 ) THEN 171: INFO = -2 172: ELSE IF( NRHS.LT.0 ) THEN 173: INFO = -3 174: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 175: INFO = -5 176: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 177: INFO = -7 178: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 179: INFO = -9 180: END IF 181: IF( INFO.NE.0 ) THEN 182: CALL XERBLA( 'DSPOSV', -INFO ) 183: RETURN 184: END IF 185: * 186: * Quick return if (N.EQ.0). 187: * 188: IF( N.EQ.0 ) 189: + RETURN 190: * 191: * Skip single precision iterative refinement if a priori slower 192: * than double precision factorization. 193: * 194: IF( .NOT.DOITREF ) THEN 195: ITER = -1 196: GO TO 40 197: END IF 198: * 199: * Compute some constants. 200: * 201: ANRM = DLANSY( 'I', UPLO, N, A, LDA, WORK ) 202: EPS = DLAMCH( 'Epsilon' ) 203: CTE = ANRM*EPS*SQRT( DBLE( N ) )*BWDMAX 204: * 205: * Set the indices PTSA, PTSX for referencing SA and SX in SWORK. 206: * 207: PTSA = 1 208: PTSX = PTSA + N*N 209: * 210: * Convert B from double precision to single precision and store the 211: * result in SX. 212: * 213: CALL DLAG2S( N, NRHS, B, LDB, SWORK( PTSX ), N, INFO ) 214: * 215: IF( INFO.NE.0 ) THEN 216: ITER = -2 217: GO TO 40 218: END IF 219: * 220: * Convert A from double precision to single precision and store the 221: * result in SA. 222: * 223: CALL DLAT2S( UPLO, N, A, LDA, SWORK( PTSA ), N, INFO ) 224: * 225: IF( INFO.NE.0 ) THEN 226: ITER = -2 227: GO TO 40 228: END IF 229: * 230: * Compute the Cholesky factorization of SA. 231: * 232: CALL SPOTRF( UPLO, N, SWORK( PTSA ), N, INFO ) 233: * 234: IF( INFO.NE.0 ) THEN 235: ITER = -3 236: GO TO 40 237: END IF 238: * 239: * Solve the system SA*SX = SB. 240: * 241: CALL SPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N, 242: + INFO ) 243: * 244: * Convert SX back to double precision 245: * 246: CALL SLAG2D( N, NRHS, SWORK( PTSX ), N, X, LDX, INFO ) 247: * 248: * Compute R = B - AX (R is WORK). 249: * 250: CALL DLACPY( 'All', N, NRHS, B, LDB, WORK, N ) 251: * 252: CALL DSYMM( 'Left', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE, 253: + WORK, N ) 254: * 255: * Check whether the NRHS normwise backward errors satisfy the 256: * stopping criterion. If yes, set ITER=0 and return. 257: * 258: DO I = 1, NRHS 259: XNRM = ABS( X( IDAMAX( N, X( 1, I ), 1 ), I ) ) 260: RNRM = ABS( WORK( IDAMAX( N, WORK( 1, I ), 1 ), I ) ) 261: IF( RNRM.GT.XNRM*CTE ) 262: + GO TO 10 263: END DO 264: * 265: * If we are here, the NRHS normwise backward errors satisfy the 266: * stopping criterion. We are good to exit. 267: * 268: ITER = 0 269: RETURN 270: * 271: 10 CONTINUE 272: * 273: DO 30 IITER = 1, ITERMAX 274: * 275: * Convert R (in WORK) from double precision to single precision 276: * and store the result in SX. 277: * 278: CALL DLAG2S( N, NRHS, WORK, N, SWORK( PTSX ), N, INFO ) 279: * 280: IF( INFO.NE.0 ) THEN 281: ITER = -2 282: GO TO 40 283: END IF 284: * 285: * Solve the system SA*SX = SR. 286: * 287: CALL SPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N, 288: + INFO ) 289: * 290: * Convert SX back to double precision and update the current 291: * iterate. 292: * 293: CALL SLAG2D( N, NRHS, SWORK( PTSX ), N, WORK, N, INFO ) 294: * 295: DO I = 1, NRHS 296: CALL DAXPY( N, ONE, WORK( 1, I ), 1, X( 1, I ), 1 ) 297: END DO 298: * 299: * Compute R = B - AX (R is WORK). 300: * 301: CALL DLACPY( 'All', N, NRHS, B, LDB, WORK, N ) 302: * 303: CALL DSYMM( 'L', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE, 304: + WORK, N ) 305: * 306: * Check whether the NRHS normwise backward errors satisfy the 307: * stopping criterion. If yes, set ITER=IITER>0 and return. 308: * 309: DO I = 1, NRHS 310: XNRM = ABS( X( IDAMAX( N, X( 1, I ), 1 ), I ) ) 311: RNRM = ABS( WORK( IDAMAX( N, WORK( 1, I ), 1 ), I ) ) 312: IF( RNRM.GT.XNRM*CTE ) 313: + GO TO 20 314: END DO 315: * 316: * If we are here, the NRHS normwise backward errors satisfy the 317: * stopping criterion, we are good to exit. 318: * 319: ITER = IITER 320: * 321: RETURN 322: * 323: 20 CONTINUE 324: * 325: 30 CONTINUE 326: * 327: * If we are at this place of the code, this is because we have 328: * performed ITER=ITERMAX iterations and never satisified the 329: * stopping criterion, set up the ITER flag accordingly and follow 330: * up on double precision routine. 331: * 332: ITER = -ITERMAX - 1 333: * 334: 40 CONTINUE 335: * 336: * Single-precision iterative refinement failed to converge to a 337: * satisfactory solution, so we resort to double precision. 338: * 339: CALL DPOTRF( UPLO, N, A, LDA, INFO ) 340: * 341: IF( INFO.NE.0 ) 342: + RETURN 343: * 344: CALL DLACPY( 'All', N, NRHS, B, LDB, X, LDX ) 345: CALL DPOTRS( UPLO, N, NRHS, A, LDA, X, LDX, INFO ) 346: * 347: RETURN 348: * 349: * End of DSPOSV. 350: * 351: END