![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DPSTRF( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO ) 2: * 3: * -- LAPACK routine (version 3.2.2) -- 4: * Craig Lucas, University of Manchester / NAG Ltd. 5: * October, 2008 6: * 7: * .. Scalar Arguments .. 8: DOUBLE PRECISION TOL 9: INTEGER INFO, LDA, N, RANK 10: CHARACTER UPLO 11: * .. 12: * .. Array Arguments .. 13: DOUBLE PRECISION A( LDA, * ), WORK( 2*N ) 14: INTEGER PIV( N ) 15: * .. 16: * 17: * Purpose 18: * ======= 19: * 20: * DPSTRF computes the Cholesky factorization with complete 21: * pivoting of a real symmetric positive semidefinite matrix A. 22: * 23: * The factorization has the form 24: * P' * A * P = U' * U , if UPLO = 'U', 25: * P' * A * P = L * L', if UPLO = 'L', 26: * where U is an upper triangular matrix and L is lower triangular, and 27: * P is stored as vector PIV. 28: * 29: * This algorithm does not attempt to check that A is positive 30: * semidefinite. This version of the algorithm calls level 3 BLAS. 31: * 32: * Arguments 33: * ========= 34: * 35: * UPLO (input) CHARACTER*1 36: * Specifies whether the upper or lower triangular part of the 37: * symmetric matrix A is stored. 38: * = 'U': Upper triangular 39: * = 'L': Lower triangular 40: * 41: * N (input) INTEGER 42: * The order of the matrix A. N >= 0. 43: * 44: * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 45: * On entry, the symmetric matrix A. If UPLO = 'U', the leading 46: * n by n upper triangular part of A contains the upper 47: * triangular part of the matrix A, and the strictly lower 48: * triangular part of A is not referenced. If UPLO = 'L', the 49: * leading n by n lower triangular part of A contains the lower 50: * triangular part of the matrix A, and the strictly upper 51: * triangular part of A is not referenced. 52: * 53: * On exit, if INFO = 0, the factor U or L from the Cholesky 54: * factorization as above. 55: * 56: * LDA (input) INTEGER 57: * The leading dimension of the array A. LDA >= max(1,N). 58: * 59: * PIV (output) INTEGER array, dimension (N) 60: * PIV is such that the nonzero entries are P( PIV(K), K ) = 1. 61: * 62: * RANK (output) INTEGER 63: * The rank of A given by the number of steps the algorithm 64: * completed. 65: * 66: * TOL (input) DOUBLE PRECISION 67: * User defined tolerance. If TOL < 0, then N*U*MAX( A(K,K) ) 68: * will be used. The algorithm terminates at the (K-1)st step 69: * if the pivot <= TOL. 70: * 71: * WORK (workspace) DOUBLE PRECISION array, dimension (2*N) 72: * Work space. 73: * 74: * INFO (output) INTEGER 75: * < 0: If INFO = -K, the K-th argument had an illegal value, 76: * = 0: algorithm completed successfully, and 77: * > 0: the matrix A is either rank deficient with computed rank 78: * as returned in RANK, or is indefinite. See Section 7 of 79: * LAPACK Working Note #161 for further information. 80: * 81: * ===================================================================== 82: * 83: * .. Parameters .. 84: DOUBLE PRECISION ONE, ZERO 85: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 86: * .. 87: * .. Local Scalars .. 88: DOUBLE PRECISION AJJ, DSTOP, DTEMP 89: INTEGER I, ITEMP, J, JB, K, NB, PVT 90: LOGICAL UPPER 91: * .. 92: * .. External Functions .. 93: DOUBLE PRECISION DLAMCH 94: INTEGER ILAENV 95: LOGICAL LSAME, DISNAN 96: EXTERNAL DLAMCH, ILAENV, LSAME, DISNAN 97: * .. 98: * .. External Subroutines .. 99: EXTERNAL DGEMV, DPSTF2, DSCAL, DSWAP, DSYRK, XERBLA 100: * .. 101: * .. Intrinsic Functions .. 102: INTRINSIC MAX, MIN, SQRT, MAXLOC 103: * .. 104: * .. Executable Statements .. 105: * 106: * Test the input parameters. 107: * 108: INFO = 0 109: UPPER = LSAME( UPLO, 'U' ) 110: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 111: INFO = -1 112: ELSE IF( N.LT.0 ) THEN 113: INFO = -2 114: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 115: INFO = -4 116: END IF 117: IF( INFO.NE.0 ) THEN 118: CALL XERBLA( 'DPSTRF', -INFO ) 119: RETURN 120: END IF 121: * 122: * Quick return if possible 123: * 124: IF( N.EQ.0 ) 125: $ RETURN 126: * 127: * Get block size 128: * 129: NB = ILAENV( 1, 'DPOTRF', UPLO, N, -1, -1, -1 ) 130: IF( NB.LE.1 .OR. NB.GE.N ) THEN 131: * 132: * Use unblocked code 133: * 134: CALL DPSTF2( UPLO, N, A( 1, 1 ), LDA, PIV, RANK, TOL, WORK, 135: $ INFO ) 136: GO TO 200 137: * 138: ELSE 139: * 140: * Initialize PIV 141: * 142: DO 100 I = 1, N 143: PIV( I ) = I 144: 100 CONTINUE 145: * 146: * Compute stopping value 147: * 148: PVT = 1 149: AJJ = A( PVT, PVT ) 150: DO I = 2, N 151: IF( A( I, I ).GT.AJJ ) THEN 152: PVT = I 153: AJJ = A( PVT, PVT ) 154: END IF 155: END DO 156: IF( AJJ.EQ.ZERO.OR.DISNAN( AJJ ) ) THEN 157: RANK = 0 158: INFO = 1 159: GO TO 200 160: END IF 161: * 162: * Compute stopping value if not supplied 163: * 164: IF( TOL.LT.ZERO ) THEN 165: DSTOP = N * DLAMCH( 'Epsilon' ) * AJJ 166: ELSE 167: DSTOP = TOL 168: END IF 169: * 170: * 171: IF( UPPER ) THEN 172: * 173: * Compute the Cholesky factorization P' * A * P = U' * U 174: * 175: DO 140 K = 1, N, NB 176: * 177: * Account for last block not being NB wide 178: * 179: JB = MIN( NB, N-K+1 ) 180: * 181: * Set relevant part of first half of WORK to zero, 182: * holds dot products 183: * 184: DO 110 I = K, N 185: WORK( I ) = 0 186: 110 CONTINUE 187: * 188: DO 130 J = K, K + JB - 1 189: * 190: * Find pivot, test for exit, else swap rows and columns 191: * Update dot products, compute possible pivots which are 192: * stored in the second half of WORK 193: * 194: DO 120 I = J, N 195: * 196: IF( J.GT.K ) THEN 197: WORK( I ) = WORK( I ) + A( J-1, I )**2 198: END IF 199: WORK( N+I ) = A( I, I ) - WORK( I ) 200: * 201: 120 CONTINUE 202: * 203: IF( J.GT.1 ) THEN 204: ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 ) 205: PVT = ITEMP + J - 1 206: AJJ = WORK( N+PVT ) 207: IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN 208: A( J, J ) = AJJ 209: GO TO 190 210: END IF 211: END IF 212: * 213: IF( J.NE.PVT ) THEN 214: * 215: * Pivot OK, so can now swap pivot rows and columns 216: * 217: A( PVT, PVT ) = A( J, J ) 218: CALL DSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 ) 219: IF( PVT.LT.N ) 220: $ CALL DSWAP( N-PVT, A( J, PVT+1 ), LDA, 221: $ A( PVT, PVT+1 ), LDA ) 222: CALL DSWAP( PVT-J-1, A( J, J+1 ), LDA, 223: $ A( J+1, PVT ), 1 ) 224: * 225: * Swap dot products and PIV 226: * 227: DTEMP = WORK( J ) 228: WORK( J ) = WORK( PVT ) 229: WORK( PVT ) = DTEMP 230: ITEMP = PIV( PVT ) 231: PIV( PVT ) = PIV( J ) 232: PIV( J ) = ITEMP 233: END IF 234: * 235: AJJ = SQRT( AJJ ) 236: A( J, J ) = AJJ 237: * 238: * Compute elements J+1:N of row J. 239: * 240: IF( J.LT.N ) THEN 241: CALL DGEMV( 'Trans', J-K, N-J, -ONE, A( K, J+1 ), 242: $ LDA, A( K, J ), 1, ONE, A( J, J+1 ), 243: $ LDA ) 244: CALL DSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA ) 245: END IF 246: * 247: 130 CONTINUE 248: * 249: * Update trailing matrix, J already incremented 250: * 251: IF( K+JB.LE.N ) THEN 252: CALL DSYRK( 'Upper', 'Trans', N-J+1, JB, -ONE, 253: $ A( K, J ), LDA, ONE, A( J, J ), LDA ) 254: END IF 255: * 256: 140 CONTINUE 257: * 258: ELSE 259: * 260: * Compute the Cholesky factorization P' * A * P = L * L' 261: * 262: DO 180 K = 1, N, NB 263: * 264: * Account for last block not being NB wide 265: * 266: JB = MIN( NB, N-K+1 ) 267: * 268: * Set relevant part of first half of WORK to zero, 269: * holds dot products 270: * 271: DO 150 I = K, N 272: WORK( I ) = 0 273: 150 CONTINUE 274: * 275: DO 170 J = K, K + JB - 1 276: * 277: * Find pivot, test for exit, else swap rows and columns 278: * Update dot products, compute possible pivots which are 279: * stored in the second half of WORK 280: * 281: DO 160 I = J, N 282: * 283: IF( J.GT.K ) THEN 284: WORK( I ) = WORK( I ) + A( I, J-1 )**2 285: END IF 286: WORK( N+I ) = A( I, I ) - WORK( I ) 287: * 288: 160 CONTINUE 289: * 290: IF( J.GT.1 ) THEN 291: ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 ) 292: PVT = ITEMP + J - 1 293: AJJ = WORK( N+PVT ) 294: IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN 295: A( J, J ) = AJJ 296: GO TO 190 297: END IF 298: END IF 299: * 300: IF( J.NE.PVT ) THEN 301: * 302: * Pivot OK, so can now swap pivot rows and columns 303: * 304: A( PVT, PVT ) = A( J, J ) 305: CALL DSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA ) 306: IF( PVT.LT.N ) 307: $ CALL DSWAP( N-PVT, A( PVT+1, J ), 1, 308: $ A( PVT+1, PVT ), 1 ) 309: CALL DSWAP( PVT-J-1, A( J+1, J ), 1, A( PVT, J+1 ), 310: $ LDA ) 311: * 312: * Swap dot products and PIV 313: * 314: DTEMP = WORK( J ) 315: WORK( J ) = WORK( PVT ) 316: WORK( PVT ) = DTEMP 317: ITEMP = PIV( PVT ) 318: PIV( PVT ) = PIV( J ) 319: PIV( J ) = ITEMP 320: END IF 321: * 322: AJJ = SQRT( AJJ ) 323: A( J, J ) = AJJ 324: * 325: * Compute elements J+1:N of column J. 326: * 327: IF( J.LT.N ) THEN 328: CALL DGEMV( 'No Trans', N-J, J-K, -ONE, 329: $ A( J+1, K ), LDA, A( J, K ), LDA, ONE, 330: $ A( J+1, J ), 1 ) 331: CALL DSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 ) 332: END IF 333: * 334: 170 CONTINUE 335: * 336: * Update trailing matrix, J already incremented 337: * 338: IF( K+JB.LE.N ) THEN 339: CALL DSYRK( 'Lower', 'No Trans', N-J+1, JB, -ONE, 340: $ A( J, K ), LDA, ONE, A( J, J ), LDA ) 341: END IF 342: * 343: 180 CONTINUE 344: * 345: END IF 346: END IF 347: * 348: * Ran to completion, A has full rank 349: * 350: RANK = N 351: * 352: GO TO 200 353: 190 CONTINUE 354: * 355: * Rank is the number of steps completed. Set INFO = 1 to signal 356: * that the factorization cannot be used to solve a system. 357: * 358: RANK = J - 1 359: INFO = 1 360: * 361: 200 CONTINUE 362: RETURN 363: * 364: * End of DPSTRF 365: * 366: END