![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DPPTRI( 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: DOUBLE PRECISION AP( * ) 14: * .. 15: * 16: * Purpose 17: * ======= 18: * 19: * DPPTRI computes the inverse of a real symmetric positive definite 20: * matrix A using the Cholesky factorization A = U**T*U or A = L*L**T 21: * computed by DPPTRF. 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) DOUBLE PRECISION array, dimension (N*(N+1)/2) 34: * On entry, the triangular factor U or L from the Cholesky 35: * factorization A = U**T*U or A = L*L**T, 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 (symmetric) 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: DOUBLE PRECISION DDOT 64: EXTERNAL LSAME, DDOT 65: * .. 66: * .. External Subroutines .. 67: EXTERNAL DSCAL, DSPR, DTPMV, DTPTRI, XERBLA 68: * .. 69: * .. Executable Statements .. 70: * 71: * Test the input parameters. 72: * 73: INFO = 0 74: UPPER = LSAME( UPLO, 'U' ) 75: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 76: INFO = -1 77: ELSE IF( N.LT.0 ) THEN 78: INFO = -2 79: END IF 80: IF( INFO.NE.0 ) THEN 81: CALL XERBLA( 'DPPTRI', -INFO ) 82: RETURN 83: END IF 84: * 85: * Quick return if possible 86: * 87: IF( N.EQ.0 ) 88: $ RETURN 89: * 90: * Invert the triangular Cholesky factor U or L. 91: * 92: CALL DTPTRI( UPLO, 'Non-unit', N, AP, INFO ) 93: IF( INFO.GT.0 ) 94: $ RETURN 95: * 96: IF( UPPER ) THEN 97: * 98: * Compute the product inv(U) * inv(U)'. 99: * 100: JJ = 0 101: DO 10 J = 1, N 102: JC = JJ + 1 103: JJ = JJ + J 104: IF( J.GT.1 ) 105: $ CALL DSPR( 'Upper', J-1, ONE, AP( JC ), 1, AP ) 106: AJJ = AP( JJ ) 107: CALL DSCAL( J, AJJ, AP( JC ), 1 ) 108: 10 CONTINUE 109: * 110: ELSE 111: * 112: * Compute the product inv(L)' * inv(L). 113: * 114: JJ = 1 115: DO 20 J = 1, N 116: JJN = JJ + N - J + 1 117: AP( JJ ) = DDOT( N-J+1, AP( JJ ), 1, AP( JJ ), 1 ) 118: IF( J.LT.N ) 119: $ CALL DTPMV( 'Lower', 'Transpose', 'Non-unit', N-J, 120: $ AP( JJN ), AP( JJ+1 ), 1 ) 121: JJ = JJN 122: 20 CONTINUE 123: END IF 124: * 125: RETURN 126: * 127: * End of DPPTRI 128: * 129: END