![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zpptri.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to zpptri.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 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