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