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