![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DPSTF2( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO ) 2: * 3: * -- LAPACK PROTOTYPE 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: * DPSTF2 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 2 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: * PIV (output) INTEGER array, dimension (N) 57: * PIV is such that the nonzero entries are P( PIV(K), K ) = 1. 58: * 59: * RANK (output) INTEGER 60: * The rank of A given by the number of steps the algorithm 61: * completed. 62: * 63: * TOL (input) DOUBLE PRECISION 64: * User defined tolerance. If TOL < 0, then N*U*MAX( A( K,K ) ) 65: * will be used. The algorithm terminates at the (K-1)st step 66: * if the pivot <= TOL. 67: * 68: * LDA (input) INTEGER 69: * The leading dimension of the array A. LDA >= max(1,N). 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, PVT 90: LOGICAL UPPER 91: * .. 92: * .. External Functions .. 93: DOUBLE PRECISION DLAMCH 94: LOGICAL LSAME, DISNAN 95: EXTERNAL DLAMCH, LSAME, DISNAN 96: * .. 97: * .. External Subroutines .. 98: EXTERNAL DGEMV, DSCAL, DSWAP, XERBLA 99: * .. 100: * .. Intrinsic Functions .. 101: INTRINSIC MAX, SQRT, MAXLOC 102: * .. 103: * .. Executable Statements .. 104: * 105: * Test the input parameters 106: * 107: INFO = 0 108: UPPER = LSAME( UPLO, 'U' ) 109: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 110: INFO = -1 111: ELSE IF( N.LT.0 ) THEN 112: INFO = -2 113: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 114: INFO = -4 115: END IF 116: IF( INFO.NE.0 ) THEN 117: CALL XERBLA( 'DPSTF2', -INFO ) 118: RETURN 119: END IF 120: * 121: * Quick return if possible 122: * 123: IF( N.EQ.0 ) 124: $ RETURN 125: * 126: * Initialize PIV 127: * 128: DO 100 I = 1, N 129: PIV( I ) = I 130: 100 CONTINUE 131: * 132: * Compute stopping value 133: * 134: PVT = 1 135: AJJ = A( PVT, PVT ) 136: DO I = 2, N 137: IF( A( I, I ).GT.AJJ ) THEN 138: PVT = I 139: AJJ = A( PVT, PVT ) 140: END IF 141: END DO 142: IF( AJJ.EQ.ZERO.OR.DISNAN( AJJ ) ) THEN 143: RANK = 0 144: INFO = 1 145: GO TO 170 146: END IF 147: * 148: * Compute stopping value if not supplied 149: * 150: IF( TOL.LT.ZERO ) THEN 151: DSTOP = N * DLAMCH( 'Epsilon' ) * AJJ 152: ELSE 153: DSTOP = TOL 154: END IF 155: * 156: * Set first half of WORK to zero, holds dot products 157: * 158: DO 110 I = 1, N 159: WORK( I ) = 0 160: 110 CONTINUE 161: * 162: IF( UPPER ) THEN 163: * 164: * Compute the Cholesky factorization P' * A * P = U' * U 165: * 166: DO 130 J = 1, N 167: * 168: * Find pivot, test for exit, else swap rows and columns 169: * Update dot products, compute possible pivots which are 170: * stored in the second half of WORK 171: * 172: DO 120 I = J, N 173: * 174: IF( J.GT.1 ) THEN 175: WORK( I ) = WORK( I ) + A( J-1, I )**2 176: END IF 177: WORK( N+I ) = A( I, I ) - WORK( I ) 178: * 179: 120 CONTINUE 180: * 181: IF( J.GT.1 ) THEN 182: ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 ) 183: PVT = ITEMP + J - 1 184: AJJ = WORK( N+PVT ) 185: IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN 186: A( J, J ) = AJJ 187: GO TO 160 188: END IF 189: END IF 190: * 191: IF( J.NE.PVT ) THEN 192: * 193: * Pivot OK, so can now swap pivot rows and columns 194: * 195: A( PVT, PVT ) = A( J, J ) 196: CALL DSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 ) 197: IF( PVT.LT.N ) 198: $ CALL DSWAP( N-PVT, A( J, PVT+1 ), LDA, 199: $ A( PVT, PVT+1 ), LDA ) 200: CALL DSWAP( PVT-J-1, A( J, J+1 ), LDA, A( J+1, PVT ), 1 ) 201: * 202: * Swap dot products and PIV 203: * 204: DTEMP = WORK( J ) 205: WORK( J ) = WORK( PVT ) 206: WORK( PVT ) = DTEMP 207: ITEMP = PIV( PVT ) 208: PIV( PVT ) = PIV( J ) 209: PIV( J ) = ITEMP 210: END IF 211: * 212: AJJ = SQRT( AJJ ) 213: A( J, J ) = AJJ 214: * 215: * Compute elements J+1:N of row J 216: * 217: IF( J.LT.N ) THEN 218: CALL DGEMV( 'Trans', J-1, N-J, -ONE, A( 1, J+1 ), LDA, 219: $ A( 1, J ), 1, ONE, A( J, J+1 ), LDA ) 220: CALL DSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA ) 221: END IF 222: * 223: 130 CONTINUE 224: * 225: ELSE 226: * 227: * Compute the Cholesky factorization P' * A * P = L * L' 228: * 229: DO 150 J = 1, N 230: * 231: * Find pivot, test for exit, else swap rows and columns 232: * Update dot products, compute possible pivots which are 233: * stored in the second half of WORK 234: * 235: DO 140 I = J, N 236: * 237: IF( J.GT.1 ) THEN 238: WORK( I ) = WORK( I ) + A( I, J-1 )**2 239: END IF 240: WORK( N+I ) = A( I, I ) - WORK( I ) 241: * 242: 140 CONTINUE 243: * 244: IF( J.GT.1 ) THEN 245: ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 ) 246: PVT = ITEMP + J - 1 247: AJJ = WORK( N+PVT ) 248: IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN 249: A( J, J ) = AJJ 250: GO TO 160 251: END IF 252: END IF 253: * 254: IF( J.NE.PVT ) THEN 255: * 256: * Pivot OK, so can now swap pivot rows and columns 257: * 258: A( PVT, PVT ) = A( J, J ) 259: CALL DSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA ) 260: IF( PVT.LT.N ) 261: $ CALL DSWAP( N-PVT, A( PVT+1, J ), 1, A( PVT+1, PVT ), 262: $ 1 ) 263: CALL DSWAP( PVT-J-1, A( J+1, J ), 1, A( PVT, J+1 ), LDA ) 264: * 265: * Swap dot products and PIV 266: * 267: DTEMP = WORK( J ) 268: WORK( J ) = WORK( PVT ) 269: WORK( PVT ) = DTEMP 270: ITEMP = PIV( PVT ) 271: PIV( PVT ) = PIV( J ) 272: PIV( J ) = ITEMP 273: END IF 274: * 275: AJJ = SQRT( AJJ ) 276: A( J, J ) = AJJ 277: * 278: * Compute elements J+1:N of column J 279: * 280: IF( J.LT.N ) THEN 281: CALL DGEMV( 'No Trans', N-J, J-1, -ONE, A( J+1, 1 ), LDA, 282: $ A( J, 1 ), LDA, ONE, A( J+1, J ), 1 ) 283: CALL DSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 ) 284: END IF 285: * 286: 150 CONTINUE 287: * 288: END IF 289: * 290: * Ran to completion, A has full rank 291: * 292: RANK = N 293: * 294: GO TO 170 295: 160 CONTINUE 296: * 297: * Rank is number of steps completed. Set INFO = 1 to signal 298: * that the factorization cannot be used to solve a system. 299: * 300: RANK = J - 1 301: INFO = 1 302: * 303: 170 CONTINUE 304: RETURN 305: * 306: * End of DPSTF2 307: * 308: END