![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZHETRF( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO ) 2: * 3: * -- LAPACK routine (version 3.2) -- 4: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 6: * November 2006 7: * 8: * .. Scalar Arguments .. 9: CHARACTER UPLO 10: INTEGER INFO, LDA, LWORK, N 11: * .. 12: * .. Array Arguments .. 13: INTEGER IPIV( * ) 14: COMPLEX*16 A( LDA, * ), WORK( * ) 15: * .. 16: * 17: * Purpose 18: * ======= 19: * 20: * ZHETRF computes the factorization of a complex Hermitian matrix A 21: * using the Bunch-Kaufman diagonal pivoting method. The form of the 22: * factorization is 23: * 24: * A = U*D*U**H or A = L*D*L**H 25: * 26: * where U (or L) is a product of permutation and unit upper (lower) 27: * triangular matrices, and D is Hermitian and block diagonal with 28: * 1-by-1 and 2-by-2 diagonal blocks. 29: * 30: * This is the blocked version of the algorithm, calling Level 3 BLAS. 31: * 32: * Arguments 33: * ========= 34: * 35: * UPLO (input) CHARACTER*1 36: * = 'U': Upper triangle of A is stored; 37: * = 'L': Lower triangle of A is stored. 38: * 39: * N (input) INTEGER 40: * The order of the matrix A. N >= 0. 41: * 42: * A (input/output) COMPLEX*16 array, dimension (LDA,N) 43: * On entry, the Hermitian matrix A. If UPLO = 'U', the leading 44: * N-by-N upper triangular part of A contains the upper 45: * triangular part of the matrix A, and the strictly lower 46: * triangular part of A is not referenced. If UPLO = 'L', the 47: * leading N-by-N lower triangular part of A contains the lower 48: * triangular part of the matrix A, and the strictly upper 49: * triangular part of A is not referenced. 50: * 51: * On exit, the block diagonal matrix D and the multipliers used 52: * to obtain the factor U or L (see below for further details). 53: * 54: * LDA (input) INTEGER 55: * The leading dimension of the array A. LDA >= max(1,N). 56: * 57: * IPIV (output) INTEGER array, dimension (N) 58: * Details of the interchanges and the block structure of D. 59: * If IPIV(k) > 0, then rows and columns k and IPIV(k) were 60: * interchanged and D(k,k) is a 1-by-1 diagonal block. 61: * If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and 62: * columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) 63: * is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = 64: * IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were 65: * interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. 66: * 67: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) 68: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 69: * 70: * LWORK (input) INTEGER 71: * The length of WORK. LWORK >=1. For best performance 72: * LWORK >= N*NB, where NB is the block size returned by ILAENV. 73: * 74: * INFO (output) INTEGER 75: * = 0: successful exit 76: * < 0: if INFO = -i, the i-th argument had an illegal value 77: * > 0: if INFO = i, D(i,i) is exactly zero. The factorization 78: * has been completed, but the block diagonal matrix D is 79: * exactly singular, and division by zero will occur if it 80: * is used to solve a system of equations. 81: * 82: * Further Details 83: * =============== 84: * 85: * If UPLO = 'U', then A = U*D*U', where 86: * U = P(n)*U(n)* ... *P(k)U(k)* ..., 87: * i.e., U is a product of terms P(k)*U(k), where k decreases from n to 88: * 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 89: * and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as 90: * defined by IPIV(k), and U(k) is a unit upper triangular matrix, such 91: * that if the diagonal block D(k) is of order s (s = 1 or 2), then 92: * 93: * ( I v 0 ) k-s 94: * U(k) = ( 0 I 0 ) s 95: * ( 0 0 I ) n-k 96: * k-s s n-k 97: * 98: * If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k). 99: * If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k), 100: * and A(k,k), and v overwrites A(1:k-2,k-1:k). 101: * 102: * If UPLO = 'L', then A = L*D*L', where 103: * L = P(1)*L(1)* ... *P(k)*L(k)* ..., 104: * i.e., L is a product of terms P(k)*L(k), where k increases from 1 to 105: * n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 106: * and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as 107: * defined by IPIV(k), and L(k) is a unit lower triangular matrix, such 108: * that if the diagonal block D(k) is of order s (s = 1 or 2), then 109: * 110: * ( I 0 0 ) k-1 111: * L(k) = ( 0 I 0 ) s 112: * ( 0 v I ) n-k-s+1 113: * k-1 s n-k-s+1 114: * 115: * If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k). 116: * If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k), 117: * and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1). 118: * 119: * ===================================================================== 120: * 121: * .. Local Scalars .. 122: LOGICAL LQUERY, UPPER 123: INTEGER IINFO, IWS, J, K, KB, LDWORK, LWKOPT, NB, NBMIN 124: * .. 125: * .. External Functions .. 126: LOGICAL LSAME 127: INTEGER ILAENV 128: EXTERNAL LSAME, ILAENV 129: * .. 130: * .. External Subroutines .. 131: EXTERNAL XERBLA, ZHETF2, ZLAHEF 132: * .. 133: * .. Intrinsic Functions .. 134: INTRINSIC MAX 135: * .. 136: * .. Executable Statements .. 137: * 138: * Test the input parameters. 139: * 140: INFO = 0 141: UPPER = LSAME( UPLO, 'U' ) 142: LQUERY = ( LWORK.EQ.-1 ) 143: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 144: INFO = -1 145: ELSE IF( N.LT.0 ) THEN 146: INFO = -2 147: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 148: INFO = -4 149: ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN 150: INFO = -7 151: END IF 152: * 153: IF( INFO.EQ.0 ) THEN 154: * 155: * Determine the block size 156: * 157: NB = ILAENV( 1, 'ZHETRF', UPLO, N, -1, -1, -1 ) 158: LWKOPT = N*NB 159: WORK( 1 ) = LWKOPT 160: END IF 161: * 162: IF( INFO.NE.0 ) THEN 163: CALL XERBLA( 'ZHETRF', -INFO ) 164: RETURN 165: ELSE IF( LQUERY ) THEN 166: RETURN 167: END IF 168: * 169: NBMIN = 2 170: LDWORK = N 171: IF( NB.GT.1 .AND. NB.LT.N ) THEN 172: IWS = LDWORK*NB 173: IF( LWORK.LT.IWS ) THEN 174: NB = MAX( LWORK / LDWORK, 1 ) 175: NBMIN = MAX( 2, ILAENV( 2, 'ZHETRF', UPLO, N, -1, -1, -1 ) ) 176: END IF 177: ELSE 178: IWS = 1 179: END IF 180: IF( NB.LT.NBMIN ) 181: $ NB = N 182: * 183: IF( UPPER ) THEN 184: * 185: * Factorize A as U*D*U' using the upper triangle of A 186: * 187: * K is the main loop index, decreasing from N to 1 in steps of 188: * KB, where KB is the number of columns factorized by ZLAHEF; 189: * KB is either NB or NB-1, or K for the last block 190: * 191: K = N 192: 10 CONTINUE 193: * 194: * If K < 1, exit from loop 195: * 196: IF( K.LT.1 ) 197: $ GO TO 40 198: * 199: IF( K.GT.NB ) THEN 200: * 201: * Factorize columns k-kb+1:k of A and use blocked code to 202: * update columns 1:k-kb 203: * 204: CALL ZLAHEF( UPLO, K, NB, KB, A, LDA, IPIV, WORK, N, IINFO ) 205: ELSE 206: * 207: * Use unblocked code to factorize columns 1:k of A 208: * 209: CALL ZHETF2( UPLO, K, A, LDA, IPIV, IINFO ) 210: KB = K 211: END IF 212: * 213: * Set INFO on the first occurrence of a zero pivot 214: * 215: IF( INFO.EQ.0 .AND. IINFO.GT.0 ) 216: $ INFO = IINFO 217: * 218: * Decrease K and return to the start of the main loop 219: * 220: K = K - KB 221: GO TO 10 222: * 223: ELSE 224: * 225: * Factorize A as L*D*L' using the lower triangle of A 226: * 227: * K is the main loop index, increasing from 1 to N in steps of 228: * KB, where KB is the number of columns factorized by ZLAHEF; 229: * KB is either NB or NB-1, or N-K+1 for the last block 230: * 231: K = 1 232: 20 CONTINUE 233: * 234: * If K > N, exit from loop 235: * 236: IF( K.GT.N ) 237: $ GO TO 40 238: * 239: IF( K.LE.N-NB ) THEN 240: * 241: * Factorize columns k:k+kb-1 of A and use blocked code to 242: * update columns k+kb:n 243: * 244: CALL ZLAHEF( UPLO, N-K+1, NB, KB, A( K, K ), LDA, IPIV( K ), 245: $ WORK, N, IINFO ) 246: ELSE 247: * 248: * Use unblocked code to factorize columns k:n of A 249: * 250: CALL ZHETF2( UPLO, N-K+1, A( K, K ), LDA, IPIV( K ), IINFO ) 251: KB = N - K + 1 252: END IF 253: * 254: * Set INFO on the first occurrence of a zero pivot 255: * 256: IF( INFO.EQ.0 .AND. IINFO.GT.0 ) 257: $ INFO = IINFO + K - 1 258: * 259: * Adjust IPIV 260: * 261: DO 30 J = K, K + KB - 1 262: IF( IPIV( J ).GT.0 ) THEN 263: IPIV( J ) = IPIV( J ) + K - 1 264: ELSE 265: IPIV( J ) = IPIV( J ) - K + 1 266: END IF 267: 30 CONTINUE 268: * 269: * Increase K and return to the start of the main loop 270: * 271: K = K + KB 272: GO TO 20 273: * 274: END IF 275: * 276: 40 CONTINUE 277: WORK( 1 ) = LWKOPT 278: RETURN 279: * 280: * End of ZHETRF 281: * 282: END