1 SUBROUTINE DPPTRI( UPLO, N, AP, INFO )
3 * -- LAPACK routine (version 3.2) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
7 * .. Scalar Arguments ..
11 * .. Array Arguments ..
12 DOUBLE PRECISION AP( * )
18 * DPPTRI computes the inverse of a real symmetric positive definite
19 * matrix A using the Cholesky factorization A = U**T*U or A = L*L**T
25 * UPLO (input) CHARACTER*1
26 * = 'U': Upper triangular factor is stored in AP;
27 * = 'L': Lower triangular factor is stored in AP.
30 * The order of the matrix A. N >= 0.
32 * AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
33 * On entry, the triangular factor U or L from the Cholesky
34 * factorization A = U**T*U or A = L*L**T, packed columnwise as
35 * a linear array. The j-th column of U or L is stored in the
36 * array AP as follows:
37 * if UPLO = 'U', AP(i + (j-1)*j/2) = U(i,j) for 1<=i<=j;
38 * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = L(i,j) for j<=i<=n.
40 * On exit, the upper or lower triangle of the (symmetric)
41 * inverse of A, overwriting the input factor U or L.
43 * INFO (output) INTEGER
44 * = 0: successful exit
45 * < 0: if INFO = -i, the i-th argument had an illegal value
46 * > 0: if INFO = i, the (i,i) element of the factor U or L is
47 * zero, and the inverse could not be computed.
49 * =====================================================================
53 PARAMETER ( ONE = 1.0D+0 )
57 INTEGER J, JC, JJ, JJN
60 * .. External Functions ..
65 * .. External Subroutines ..
66 EXTERNAL DSCAL, DSPR, DTPMV, DTPTRI, XERBLA
68 * .. Executable Statements ..
70 * Test the input parameters.
73 UPPER = LSAME( UPLO, 'U' )
74 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
76 ELSE IF( N.LT.0 ) THEN
80 CALL XERBLA( 'DPPTRI', -INFO )
84 * Quick return if possible
89 * Invert the triangular Cholesky factor U or L.
91 CALL DTPTRI( UPLO, 'Non-unit', N, AP, INFO )
97 * Compute the product inv(U) * inv(U)'.
104 $ CALL DSPR( 'Upper', J-1, ONE, AP( JC ), 1, AP )
106 CALL DSCAL( J, AJJ, AP( JC ), 1 )
111 * Compute the product inv(L)' * inv(L).
116 AP( JJ ) = DDOT( N-J+1, AP( JJ ), 1, AP( JJ ), 1 )
118 $ CALL DTPMV( 'Lower', 'Transpose', 'Non-unit', N-J,
119 $ AP( JJN ), AP( JJ+1 ), 1 )