![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZHETRI( UPLO, N, A, LDA, IPIV, WORK, 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, N 11: * .. 12: * .. Array Arguments .. 13: INTEGER IPIV( * ) 14: COMPLEX*16 A( LDA, * ), WORK( * ) 15: * .. 16: * 17: * Purpose 18: * ======= 19: * 20: * ZHETRI computes the inverse of a complex Hermitian indefinite matrix 21: * A using the factorization A = U*D*U**H or A = L*D*L**H computed by 22: * ZHETRF. 23: * 24: * Arguments 25: * ========= 26: * 27: * UPLO (input) CHARACTER*1 28: * Specifies whether the details of the factorization are stored 29: * as an upper or lower triangular matrix. 30: * = 'U': Upper triangular, form is A = U*D*U**H; 31: * = 'L': Lower triangular, form is A = L*D*L**H. 32: * 33: * N (input) INTEGER 34: * The order of the matrix A. N >= 0. 35: * 36: * A (input/output) COMPLEX*16 array, dimension (LDA,N) 37: * On entry, the block diagonal matrix D and the multipliers 38: * used to obtain the factor U or L as computed by ZHETRF. 39: * 40: * On exit, if INFO = 0, the (Hermitian) inverse of the original 41: * matrix. If UPLO = 'U', the upper triangular part of the 42: * inverse is formed and the part of A below the diagonal is not 43: * referenced; if UPLO = 'L' the lower triangular part of the 44: * inverse is formed and the part of A above the diagonal is 45: * not referenced. 46: * 47: * LDA (input) INTEGER 48: * The leading dimension of the array A. LDA >= max(1,N). 49: * 50: * IPIV (input) INTEGER array, dimension (N) 51: * Details of the interchanges and the block structure of D 52: * as determined by ZHETRF. 53: * 54: * WORK (workspace) COMPLEX*16 array, dimension (N) 55: * 56: * INFO (output) INTEGER 57: * = 0: successful exit 58: * < 0: if INFO = -i, the i-th argument had an illegal value 59: * > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its 60: * inverse could not be computed. 61: * 62: * ===================================================================== 63: * 64: * .. Parameters .. 65: DOUBLE PRECISION ONE 66: COMPLEX*16 CONE, ZERO 67: PARAMETER ( ONE = 1.0D+0, CONE = ( 1.0D+0, 0.0D+0 ), 68: $ ZERO = ( 0.0D+0, 0.0D+0 ) ) 69: * .. 70: * .. Local Scalars .. 71: LOGICAL UPPER 72: INTEGER J, K, KP, KSTEP 73: DOUBLE PRECISION AK, AKP1, D, T 74: COMPLEX*16 AKKP1, TEMP 75: * .. 76: * .. External Functions .. 77: LOGICAL LSAME 78: COMPLEX*16 ZDOTC 79: EXTERNAL LSAME, ZDOTC 80: * .. 81: * .. External Subroutines .. 82: EXTERNAL XERBLA, ZCOPY, ZHEMV, ZSWAP 83: * .. 84: * .. Intrinsic Functions .. 85: INTRINSIC ABS, DBLE, DCONJG, MAX 86: * .. 87: * .. Executable Statements .. 88: * 89: * Test the input parameters. 90: * 91: INFO = 0 92: UPPER = LSAME( UPLO, 'U' ) 93: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 94: INFO = -1 95: ELSE IF( N.LT.0 ) THEN 96: INFO = -2 97: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 98: INFO = -4 99: END IF 100: IF( INFO.NE.0 ) THEN 101: CALL XERBLA( 'ZHETRI', -INFO ) 102: RETURN 103: END IF 104: * 105: * Quick return if possible 106: * 107: IF( N.EQ.0 ) 108: $ RETURN 109: * 110: * Check that the diagonal matrix D is nonsingular. 111: * 112: IF( UPPER ) THEN 113: * 114: * Upper triangular storage: examine D from bottom to top 115: * 116: DO 10 INFO = N, 1, -1 117: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO ) 118: $ RETURN 119: 10 CONTINUE 120: ELSE 121: * 122: * Lower triangular storage: examine D from top to bottom. 123: * 124: DO 20 INFO = 1, N 125: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO ) 126: $ RETURN 127: 20 CONTINUE 128: END IF 129: INFO = 0 130: * 131: IF( UPPER ) THEN 132: * 133: * Compute inv(A) from the factorization A = U*D*U'. 134: * 135: * K is the main loop index, increasing from 1 to N in steps of 136: * 1 or 2, depending on the size of the diagonal blocks. 137: * 138: K = 1 139: 30 CONTINUE 140: * 141: * If K > N, exit from loop. 142: * 143: IF( K.GT.N ) 144: $ GO TO 50 145: * 146: IF( IPIV( K ).GT.0 ) THEN 147: * 148: * 1 x 1 diagonal block 149: * 150: * Invert the diagonal block. 151: * 152: A( K, K ) = ONE / DBLE( A( K, K ) ) 153: * 154: * Compute column K of the inverse. 155: * 156: IF( K.GT.1 ) THEN 157: CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 ) 158: CALL ZHEMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, ZERO, 159: $ A( 1, K ), 1 ) 160: A( K, K ) = A( K, K ) - DBLE( ZDOTC( K-1, WORK, 1, A( 1, 161: $ K ), 1 ) ) 162: END IF 163: KSTEP = 1 164: ELSE 165: * 166: * 2 x 2 diagonal block 167: * 168: * Invert the diagonal block. 169: * 170: T = ABS( A( K, K+1 ) ) 171: AK = DBLE( A( K, K ) ) / T 172: AKP1 = DBLE( A( K+1, K+1 ) ) / T 173: AKKP1 = A( K, K+1 ) / T 174: D = T*( AK*AKP1-ONE ) 175: A( K, K ) = AKP1 / D 176: A( K+1, K+1 ) = AK / D 177: A( K, K+1 ) = -AKKP1 / D 178: * 179: * Compute columns K and K+1 of the inverse. 180: * 181: IF( K.GT.1 ) THEN 182: CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 ) 183: CALL ZHEMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, ZERO, 184: $ A( 1, K ), 1 ) 185: A( K, K ) = A( K, K ) - DBLE( ZDOTC( K-1, WORK, 1, A( 1, 186: $ K ), 1 ) ) 187: A( K, K+1 ) = A( K, K+1 ) - 188: $ ZDOTC( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 ) 189: CALL ZCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 ) 190: CALL ZHEMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, ZERO, 191: $ A( 1, K+1 ), 1 ) 192: A( K+1, K+1 ) = A( K+1, K+1 ) - 193: $ DBLE( ZDOTC( K-1, WORK, 1, A( 1, K+1 ), 194: $ 1 ) ) 195: END IF 196: KSTEP = 2 197: END IF 198: * 199: KP = ABS( IPIV( K ) ) 200: IF( KP.NE.K ) THEN 201: * 202: * Interchange rows and columns K and KP in the leading 203: * submatrix A(1:k+1,1:k+1) 204: * 205: CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 ) 206: DO 40 J = KP + 1, K - 1 207: TEMP = DCONJG( A( J, K ) ) 208: A( J, K ) = DCONJG( A( KP, J ) ) 209: A( KP, J ) = TEMP 210: 40 CONTINUE 211: A( KP, K ) = DCONJG( A( KP, K ) ) 212: TEMP = A( K, K ) 213: A( K, K ) = A( KP, KP ) 214: A( KP, KP ) = TEMP 215: IF( KSTEP.EQ.2 ) THEN 216: TEMP = A( K, K+1 ) 217: A( K, K+1 ) = A( KP, K+1 ) 218: A( KP, K+1 ) = TEMP 219: END IF 220: END IF 221: * 222: K = K + KSTEP 223: GO TO 30 224: 50 CONTINUE 225: * 226: ELSE 227: * 228: * Compute inv(A) from the factorization A = L*D*L'. 229: * 230: * K is the main loop index, increasing from 1 to N in steps of 231: * 1 or 2, depending on the size of the diagonal blocks. 232: * 233: K = N 234: 60 CONTINUE 235: * 236: * If K < 1, exit from loop. 237: * 238: IF( K.LT.1 ) 239: $ GO TO 80 240: * 241: IF( IPIV( K ).GT.0 ) THEN 242: * 243: * 1 x 1 diagonal block 244: * 245: * Invert the diagonal block. 246: * 247: A( K, K ) = ONE / DBLE( A( K, K ) ) 248: * 249: * Compute column K of the inverse. 250: * 251: IF( K.LT.N ) THEN 252: CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 ) 253: CALL ZHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK, 254: $ 1, ZERO, A( K+1, K ), 1 ) 255: A( K, K ) = A( K, K ) - DBLE( ZDOTC( N-K, WORK, 1, 256: $ A( K+1, K ), 1 ) ) 257: END IF 258: KSTEP = 1 259: ELSE 260: * 261: * 2 x 2 diagonal block 262: * 263: * Invert the diagonal block. 264: * 265: T = ABS( A( K, K-1 ) ) 266: AK = DBLE( A( K-1, K-1 ) ) / T 267: AKP1 = DBLE( A( K, K ) ) / T 268: AKKP1 = A( K, K-1 ) / T 269: D = T*( AK*AKP1-ONE ) 270: A( K-1, K-1 ) = AKP1 / D 271: A( K, K ) = AK / D 272: A( K, K-1 ) = -AKKP1 / D 273: * 274: * Compute columns K-1 and K of the inverse. 275: * 276: IF( K.LT.N ) THEN 277: CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 ) 278: CALL ZHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK, 279: $ 1, ZERO, A( K+1, K ), 1 ) 280: A( K, K ) = A( K, K ) - DBLE( ZDOTC( N-K, WORK, 1, 281: $ A( K+1, K ), 1 ) ) 282: A( K, K-1 ) = A( K, K-1 ) - 283: $ ZDOTC( N-K, A( K+1, K ), 1, A( K+1, K-1 ), 284: $ 1 ) 285: CALL ZCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 ) 286: CALL ZHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK, 287: $ 1, ZERO, A( K+1, K-1 ), 1 ) 288: A( K-1, K-1 ) = A( K-1, K-1 ) - 289: $ DBLE( ZDOTC( N-K, WORK, 1, A( K+1, K-1 ), 290: $ 1 ) ) 291: END IF 292: KSTEP = 2 293: END IF 294: * 295: KP = ABS( IPIV( K ) ) 296: IF( KP.NE.K ) THEN 297: * 298: * Interchange rows and columns K and KP in the trailing 299: * submatrix A(k-1:n,k-1:n) 300: * 301: IF( KP.LT.N ) 302: $ CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 ) 303: DO 70 J = K + 1, KP - 1 304: TEMP = DCONJG( A( J, K ) ) 305: A( J, K ) = DCONJG( A( KP, J ) ) 306: A( KP, J ) = TEMP 307: 70 CONTINUE 308: A( KP, K ) = DCONJG( A( KP, K ) ) 309: TEMP = A( K, K ) 310: A( K, K ) = A( KP, KP ) 311: A( KP, KP ) = TEMP 312: IF( KSTEP.EQ.2 ) THEN 313: TEMP = A( K, K-1 ) 314: A( K, K-1 ) = A( KP, K-1 ) 315: A( KP, K-1 ) = TEMP 316: END IF 317: END IF 318: * 319: K = K - KSTEP 320: GO TO 60 321: 80 CONTINUE 322: END IF 323: * 324: RETURN 325: * 326: * End of ZHETRI 327: * 328: END