![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZPPTRI( UPLO, N, AP, 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: COMPLEX*16 AP( * ) 14: * .. 15: * 16: * Purpose 17: * ======= 18: * 19: * ZPPTRI computes the inverse of a complex Hermitian positive definite 20: * matrix A using the Cholesky factorization A = U**H*U or A = L*L**H 21: * computed by ZPPTRF. 22: * 23: * Arguments 24: * ========= 25: * 26: * UPLO (input) CHARACTER*1 27: * = 'U': Upper triangular factor is stored in AP; 28: * = 'L': Lower triangular factor is stored in AP. 29: * 30: * N (input) INTEGER 31: * The order of the matrix A. N >= 0. 32: * 33: * AP (input/output) COMPLEX*16 array, dimension (N*(N+1)/2) 34: * On entry, the triangular factor U or L from the Cholesky 35: * factorization A = U**H*U or A = L*L**H, packed columnwise as 36: * a linear array. The j-th column of U or L is stored in the 37: * array AP as follows: 38: * if UPLO = 'U', AP(i + (j-1)*j/2) = U(i,j) for 1<=i<=j; 39: * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = L(i,j) for j<=i<=n. 40: * 41: * On exit, the upper or lower triangle of the (Hermitian) 42: * inverse of A, overwriting the input factor U or L. 43: * 44: * INFO (output) INTEGER 45: * = 0: successful exit 46: * < 0: if INFO = -i, the i-th argument had an illegal value 47: * > 0: if INFO = i, the (i,i) element of the factor U or L is 48: * zero, and the inverse could not be computed. 49: * 50: * ===================================================================== 51: * 52: * .. Parameters .. 53: DOUBLE PRECISION ONE 54: PARAMETER ( ONE = 1.0D+0 ) 55: * .. 56: * .. Local Scalars .. 57: LOGICAL UPPER 58: INTEGER J, JC, JJ, JJN 59: DOUBLE PRECISION AJJ 60: * .. 61: * .. External Functions .. 62: LOGICAL LSAME 63: COMPLEX*16 ZDOTC 64: EXTERNAL LSAME, ZDOTC 65: * .. 66: * .. External Subroutines .. 67: EXTERNAL XERBLA, ZDSCAL, ZHPR, ZTPMV, ZTPTRI 68: * .. 69: * .. Intrinsic Functions .. 70: INTRINSIC DBLE 71: * .. 72: * .. Executable Statements .. 73: * 74: * Test the input parameters. 75: * 76: INFO = 0 77: UPPER = LSAME( UPLO, 'U' ) 78: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 79: INFO = -1 80: ELSE IF( N.LT.0 ) THEN 81: INFO = -2 82: END IF 83: IF( INFO.NE.0 ) THEN 84: CALL XERBLA( 'ZPPTRI', -INFO ) 85: RETURN 86: END IF 87: * 88: * Quick return if possible 89: * 90: IF( N.EQ.0 ) 91: $ RETURN 92: * 93: * Invert the triangular Cholesky factor U or L. 94: * 95: CALL ZTPTRI( UPLO, 'Non-unit', N, AP, INFO ) 96: IF( INFO.GT.0 ) 97: $ RETURN 98: IF( UPPER ) THEN 99: * 100: * Compute the product inv(U) * inv(U)'. 101: * 102: JJ = 0 103: DO 10 J = 1, N 104: JC = JJ + 1 105: JJ = JJ + J 106: IF( J.GT.1 ) 107: $ CALL ZHPR( 'Upper', J-1, ONE, AP( JC ), 1, AP ) 108: AJJ = AP( JJ ) 109: CALL ZDSCAL( J, AJJ, AP( JC ), 1 ) 110: 10 CONTINUE 111: * 112: ELSE 113: * 114: * Compute the product inv(L)' * inv(L). 115: * 116: JJ = 1 117: DO 20 J = 1, N 118: JJN = JJ + N - J + 1 119: AP( JJ ) = DBLE( ZDOTC( N-J+1, AP( JJ ), 1, AP( JJ ), 1 ) ) 120: IF( J.LT.N ) 121: $ CALL ZTPMV( 'Lower', 'Conjugate transpose', 'Non-unit', 122: $ N-J, AP( JJN ), AP( JJ+1 ), 1 ) 123: JJ = JJN 124: 20 CONTINUE 125: END IF 126: * 127: RETURN 128: * 129: * End of ZPPTRI 130: * 131: END