![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZSYTRI( 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: * ZSYTRI computes the inverse of a complex symmetric indefinite matrix 21: * A using the factorization A = U*D*U**T or A = L*D*L**T computed by 22: * ZSYTRF. 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**T; 31: * = 'L': Lower triangular, form is A = L*D*L**T. 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 ZSYTRF. 39: * 40: * On exit, if INFO = 0, the (symmetric) 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 ZSYTRF. 53: * 54: * WORK (workspace) COMPLEX*16 array, dimension (2*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: COMPLEX*16 ONE, ZERO 66: PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ), 67: $ ZERO = ( 0.0D+0, 0.0D+0 ) ) 68: * .. 69: * .. Local Scalars .. 70: LOGICAL UPPER 71: INTEGER K, KP, KSTEP 72: COMPLEX*16 AK, AKKP1, AKP1, D, T, TEMP 73: * .. 74: * .. External Functions .. 75: LOGICAL LSAME 76: COMPLEX*16 ZDOTU 77: EXTERNAL LSAME, ZDOTU 78: * .. 79: * .. External Subroutines .. 80: EXTERNAL XERBLA, ZCOPY, ZSWAP, ZSYMV 81: * .. 82: * .. Intrinsic Functions .. 83: INTRINSIC ABS, MAX 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: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 96: INFO = -4 97: END IF 98: IF( INFO.NE.0 ) THEN 99: CALL XERBLA( 'ZSYTRI', -INFO ) 100: RETURN 101: END IF 102: * 103: * Quick return if possible 104: * 105: IF( N.EQ.0 ) 106: $ RETURN 107: * 108: * Check that the diagonal matrix D is nonsingular. 109: * 110: IF( UPPER ) THEN 111: * 112: * Upper triangular storage: examine D from bottom to top 113: * 114: DO 10 INFO = N, 1, -1 115: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO ) 116: $ RETURN 117: 10 CONTINUE 118: ELSE 119: * 120: * Lower triangular storage: examine D from top to bottom. 121: * 122: DO 20 INFO = 1, N 123: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO ) 124: $ RETURN 125: 20 CONTINUE 126: END IF 127: INFO = 0 128: * 129: IF( UPPER ) THEN 130: * 131: * Compute inv(A) from the factorization A = U*D*U'. 132: * 133: * K is the main loop index, increasing from 1 to N in steps of 134: * 1 or 2, depending on the size of the diagonal blocks. 135: * 136: K = 1 137: 30 CONTINUE 138: * 139: * If K > N, exit from loop. 140: * 141: IF( K.GT.N ) 142: $ GO TO 40 143: * 144: IF( IPIV( K ).GT.0 ) THEN 145: * 146: * 1 x 1 diagonal block 147: * 148: * Invert the diagonal block. 149: * 150: A( K, K ) = ONE / A( K, K ) 151: * 152: * Compute column K of the inverse. 153: * 154: IF( K.GT.1 ) THEN 155: CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 ) 156: CALL ZSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO, 157: $ A( 1, K ), 1 ) 158: A( K, K ) = A( K, K ) - ZDOTU( K-1, WORK, 1, A( 1, K ), 159: $ 1 ) 160: END IF 161: KSTEP = 1 162: ELSE 163: * 164: * 2 x 2 diagonal block 165: * 166: * Invert the diagonal block. 167: * 168: T = A( K, K+1 ) 169: AK = A( K, K ) / T 170: AKP1 = A( K+1, K+1 ) / T 171: AKKP1 = A( K, K+1 ) / T 172: D = T*( AK*AKP1-ONE ) 173: A( K, K ) = AKP1 / D 174: A( K+1, K+1 ) = AK / D 175: A( K, K+1 ) = -AKKP1 / D 176: * 177: * Compute columns K and K+1 of the inverse. 178: * 179: IF( K.GT.1 ) THEN 180: CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 ) 181: CALL ZSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO, 182: $ A( 1, K ), 1 ) 183: A( K, K ) = A( K, K ) - ZDOTU( K-1, WORK, 1, A( 1, K ), 184: $ 1 ) 185: A( K, K+1 ) = A( K, K+1 ) - 186: $ ZDOTU( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 ) 187: CALL ZCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 ) 188: CALL ZSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO, 189: $ A( 1, K+1 ), 1 ) 190: A( K+1, K+1 ) = A( K+1, K+1 ) - 191: $ ZDOTU( K-1, WORK, 1, A( 1, K+1 ), 1 ) 192: END IF 193: KSTEP = 2 194: END IF 195: * 196: KP = ABS( IPIV( K ) ) 197: IF( KP.NE.K ) THEN 198: * 199: * Interchange rows and columns K and KP in the leading 200: * submatrix A(1:k+1,1:k+1) 201: * 202: CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 ) 203: CALL ZSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA ) 204: TEMP = A( K, K ) 205: A( K, K ) = A( KP, KP ) 206: A( KP, KP ) = TEMP 207: IF( KSTEP.EQ.2 ) THEN 208: TEMP = A( K, K+1 ) 209: A( K, K+1 ) = A( KP, K+1 ) 210: A( KP, K+1 ) = TEMP 211: END IF 212: END IF 213: * 214: K = K + KSTEP 215: GO TO 30 216: 40 CONTINUE 217: * 218: ELSE 219: * 220: * Compute inv(A) from the factorization A = L*D*L'. 221: * 222: * K is the main loop index, increasing from 1 to N in steps of 223: * 1 or 2, depending on the size of the diagonal blocks. 224: * 225: K = N 226: 50 CONTINUE 227: * 228: * If K < 1, exit from loop. 229: * 230: IF( K.LT.1 ) 231: $ GO TO 60 232: * 233: IF( IPIV( K ).GT.0 ) THEN 234: * 235: * 1 x 1 diagonal block 236: * 237: * Invert the diagonal block. 238: * 239: A( K, K ) = ONE / A( K, K ) 240: * 241: * Compute column K of the inverse. 242: * 243: IF( K.LT.N ) THEN 244: CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 ) 245: CALL ZSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1, 246: $ ZERO, A( K+1, K ), 1 ) 247: A( K, K ) = A( K, K ) - ZDOTU( N-K, WORK, 1, A( K+1, K ), 248: $ 1 ) 249: END IF 250: KSTEP = 1 251: ELSE 252: * 253: * 2 x 2 diagonal block 254: * 255: * Invert the diagonal block. 256: * 257: T = A( K, K-1 ) 258: AK = A( K-1, K-1 ) / T 259: AKP1 = A( K, K ) / T 260: AKKP1 = A( K, K-1 ) / T 261: D = T*( AK*AKP1-ONE ) 262: A( K-1, K-1 ) = AKP1 / D 263: A( K, K ) = AK / D 264: A( K, K-1 ) = -AKKP1 / D 265: * 266: * Compute columns K-1 and K of the inverse. 267: * 268: IF( K.LT.N ) THEN 269: CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 ) 270: CALL ZSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1, 271: $ ZERO, A( K+1, K ), 1 ) 272: A( K, K ) = A( K, K ) - ZDOTU( N-K, WORK, 1, A( K+1, K ), 273: $ 1 ) 274: A( K, K-1 ) = A( K, K-1 ) - 275: $ ZDOTU( N-K, A( K+1, K ), 1, A( K+1, K-1 ), 276: $ 1 ) 277: CALL ZCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 ) 278: CALL ZSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1, 279: $ ZERO, A( K+1, K-1 ), 1 ) 280: A( K-1, K-1 ) = A( K-1, K-1 ) - 281: $ ZDOTU( N-K, WORK, 1, A( K+1, K-1 ), 1 ) 282: END IF 283: KSTEP = 2 284: END IF 285: * 286: KP = ABS( IPIV( K ) ) 287: IF( KP.NE.K ) THEN 288: * 289: * Interchange rows and columns K and KP in the trailing 290: * submatrix A(k-1:n,k-1:n) 291: * 292: IF( KP.LT.N ) 293: $ CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 ) 294: CALL ZSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA ) 295: TEMP = A( K, K ) 296: A( K, K ) = A( KP, KP ) 297: A( KP, KP ) = TEMP 298: IF( KSTEP.EQ.2 ) THEN 299: TEMP = A( K, K-1 ) 300: A( K, K-1 ) = A( KP, K-1 ) 301: A( KP, K-1 ) = TEMP 302: END IF 303: END IF 304: * 305: K = K - KSTEP 306: GO TO 50 307: 60 CONTINUE 308: END IF 309: * 310: RETURN 311: * 312: * End of ZSYTRI 313: * 314: END