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
CVSweb interface <joel.bertrand@systella.fr>