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