1 SUBROUTINE ZPPTRF( UPLO, N, AP, INFO )
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
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 * ZPPTRF computes the Cholesky factorization of a complex Hermitian
20 * positive definite matrix A stored in packed format.
21 *
22 * The factorization has the form
23 * A = U**H * U, if UPLO = 'U', or
24 * A = L * L**H, if UPLO = 'L',
25 * where U is an upper triangular matrix and L is lower triangular.
26 *
27 * Arguments
28 * =========
29 *
30 * UPLO (input) CHARACTER*1
31 * = 'U': Upper triangle of A is stored;
32 * = 'L': Lower triangle of A is stored.
33 *
34 * N (input) INTEGER
35 * The order of the matrix A. N >= 0.
36 *
37 * AP (input/output) COMPLEX*16 array, dimension (N*(N+1)/2)
38 * On entry, the upper or lower triangle of the Hermitian matrix
39 * A, packed columnwise in a linear array. The j-th column of A
40 * is stored in the array AP as follows:
41 * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
42 * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
43 * See below for further details.
44 *
45 * On exit, if INFO = 0, the triangular factor U or L from the
46 * Cholesky factorization A = U**H*U or A = L*L**H, in the same
47 * storage format as A.
48 *
49 * INFO (output) INTEGER
50 * = 0: successful exit
51 * < 0: if INFO = -i, the i-th argument had an illegal value
52 * > 0: if INFO = i, the leading minor of order i is not
53 * positive definite, and the factorization could not be
54 * completed.
55 *
56 * Further Details
57 * ===============
58 *
59 * The packed storage scheme is illustrated by the following example
60 * when N = 4, UPLO = 'U':
61 *
62 * Two-dimensional storage of the Hermitian matrix A:
63 *
64 * a11 a12 a13 a14
65 * a22 a23 a24
66 * a33 a34 (aij = conjg(aji))
67 * a44
68 *
69 * Packed storage of the upper triangle of A:
70 *
71 * AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
72 *
73 * =====================================================================
74 *
75 * .. Parameters ..
76 DOUBLE PRECISION ZERO, ONE
77 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
78 * ..
79 * .. Local Scalars ..
80 LOGICAL UPPER
81 INTEGER J, JC, JJ
82 DOUBLE PRECISION AJJ
83 * ..
84 * .. External Functions ..
85 LOGICAL LSAME
86 COMPLEX*16 ZDOTC
87 EXTERNAL LSAME, ZDOTC
88 * ..
89 * .. External Subroutines ..
90 EXTERNAL XERBLA, ZDSCAL, ZHPR, ZTPSV
91 * ..
92 * .. Intrinsic Functions ..
93 INTRINSIC DBLE, SQRT
94 * ..
95 * .. Executable Statements ..
96 *
97 * Test the input parameters.
98 *
99 INFO = 0
100 UPPER = LSAME( UPLO, 'U' )
101 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
102 INFO = -1
103 ELSE IF( N.LT.0 ) THEN
104 INFO = -2
105 END IF
106 IF( INFO.NE.0 ) THEN
107 CALL XERBLA( 'ZPPTRF', -INFO )
108 RETURN
109 END IF
110 *
111 * Quick return if possible
112 *
113 IF( N.EQ.0 )
114 $ RETURN
115 *
116 IF( UPPER ) THEN
117 *
118 * Compute the Cholesky factorization A = U**H * U.
119 *
120 JJ = 0
121 DO 10 J = 1, N
122 JC = JJ + 1
123 JJ = JJ + J
124 *
125 * Compute elements 1:J-1 of column J.
126 *
127 IF( J.GT.1 )
128 $ CALL ZTPSV( 'Upper', 'Conjugate transpose', 'Non-unit',
129 $ J-1, AP, AP( JC ), 1 )
130 *
131 * Compute U(J,J) and test for non-positive-definiteness.
132 *
133 AJJ = DBLE( AP( JJ ) ) - ZDOTC( J-1, AP( JC ), 1, AP( JC ),
134 $ 1 )
135 IF( AJJ.LE.ZERO ) THEN
136 AP( JJ ) = AJJ
137 GO TO 30
138 END IF
139 AP( JJ ) = SQRT( AJJ )
140 10 CONTINUE
141 ELSE
142 *
143 * Compute the Cholesky factorization A = L * L**H.
144 *
145 JJ = 1
146 DO 20 J = 1, N
147 *
148 * Compute L(J,J) and test for non-positive-definiteness.
149 *
150 AJJ = DBLE( AP( JJ ) )
151 IF( AJJ.LE.ZERO ) THEN
152 AP( JJ ) = AJJ
153 GO TO 30
154 END IF
155 AJJ = SQRT( AJJ )
156 AP( JJ ) = AJJ
157 *
158 * Compute elements J+1:N of column J and update the trailing
159 * submatrix.
160 *
161 IF( J.LT.N ) THEN
162 CALL ZDSCAL( N-J, ONE / AJJ, AP( JJ+1 ), 1 )
163 CALL ZHPR( 'Lower', N-J, -ONE, AP( JJ+1 ), 1,
164 $ AP( JJ+N-J+1 ) )
165 JJ = JJ + N - J + 1
166 END IF
167 20 CONTINUE
168 END IF
169 GO TO 40
170 *
171 30 CONTINUE
172 INFO = J
173 *
174 40 CONTINUE
175 RETURN
176 *
177 * End of ZPPTRF
178 *
179 END
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
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 * ZPPTRF computes the Cholesky factorization of a complex Hermitian
20 * positive definite matrix A stored in packed format.
21 *
22 * The factorization has the form
23 * A = U**H * U, if UPLO = 'U', or
24 * A = L * L**H, if UPLO = 'L',
25 * where U is an upper triangular matrix and L is lower triangular.
26 *
27 * Arguments
28 * =========
29 *
30 * UPLO (input) CHARACTER*1
31 * = 'U': Upper triangle of A is stored;
32 * = 'L': Lower triangle of A is stored.
33 *
34 * N (input) INTEGER
35 * The order of the matrix A. N >= 0.
36 *
37 * AP (input/output) COMPLEX*16 array, dimension (N*(N+1)/2)
38 * On entry, the upper or lower triangle of the Hermitian matrix
39 * A, packed columnwise in a linear array. The j-th column of A
40 * is stored in the array AP as follows:
41 * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
42 * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
43 * See below for further details.
44 *
45 * On exit, if INFO = 0, the triangular factor U or L from the
46 * Cholesky factorization A = U**H*U or A = L*L**H, in the same
47 * storage format as A.
48 *
49 * INFO (output) INTEGER
50 * = 0: successful exit
51 * < 0: if INFO = -i, the i-th argument had an illegal value
52 * > 0: if INFO = i, the leading minor of order i is not
53 * positive definite, and the factorization could not be
54 * completed.
55 *
56 * Further Details
57 * ===============
58 *
59 * The packed storage scheme is illustrated by the following example
60 * when N = 4, UPLO = 'U':
61 *
62 * Two-dimensional storage of the Hermitian matrix A:
63 *
64 * a11 a12 a13 a14
65 * a22 a23 a24
66 * a33 a34 (aij = conjg(aji))
67 * a44
68 *
69 * Packed storage of the upper triangle of A:
70 *
71 * AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
72 *
73 * =====================================================================
74 *
75 * .. Parameters ..
76 DOUBLE PRECISION ZERO, ONE
77 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
78 * ..
79 * .. Local Scalars ..
80 LOGICAL UPPER
81 INTEGER J, JC, JJ
82 DOUBLE PRECISION AJJ
83 * ..
84 * .. External Functions ..
85 LOGICAL LSAME
86 COMPLEX*16 ZDOTC
87 EXTERNAL LSAME, ZDOTC
88 * ..
89 * .. External Subroutines ..
90 EXTERNAL XERBLA, ZDSCAL, ZHPR, ZTPSV
91 * ..
92 * .. Intrinsic Functions ..
93 INTRINSIC DBLE, SQRT
94 * ..
95 * .. Executable Statements ..
96 *
97 * Test the input parameters.
98 *
99 INFO = 0
100 UPPER = LSAME( UPLO, 'U' )
101 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
102 INFO = -1
103 ELSE IF( N.LT.0 ) THEN
104 INFO = -2
105 END IF
106 IF( INFO.NE.0 ) THEN
107 CALL XERBLA( 'ZPPTRF', -INFO )
108 RETURN
109 END IF
110 *
111 * Quick return if possible
112 *
113 IF( N.EQ.0 )
114 $ RETURN
115 *
116 IF( UPPER ) THEN
117 *
118 * Compute the Cholesky factorization A = U**H * U.
119 *
120 JJ = 0
121 DO 10 J = 1, N
122 JC = JJ + 1
123 JJ = JJ + J
124 *
125 * Compute elements 1:J-1 of column J.
126 *
127 IF( J.GT.1 )
128 $ CALL ZTPSV( 'Upper', 'Conjugate transpose', 'Non-unit',
129 $ J-1, AP, AP( JC ), 1 )
130 *
131 * Compute U(J,J) and test for non-positive-definiteness.
132 *
133 AJJ = DBLE( AP( JJ ) ) - ZDOTC( J-1, AP( JC ), 1, AP( JC ),
134 $ 1 )
135 IF( AJJ.LE.ZERO ) THEN
136 AP( JJ ) = AJJ
137 GO TO 30
138 END IF
139 AP( JJ ) = SQRT( AJJ )
140 10 CONTINUE
141 ELSE
142 *
143 * Compute the Cholesky factorization A = L * L**H.
144 *
145 JJ = 1
146 DO 20 J = 1, N
147 *
148 * Compute L(J,J) and test for non-positive-definiteness.
149 *
150 AJJ = DBLE( AP( JJ ) )
151 IF( AJJ.LE.ZERO ) THEN
152 AP( JJ ) = AJJ
153 GO TO 30
154 END IF
155 AJJ = SQRT( AJJ )
156 AP( JJ ) = AJJ
157 *
158 * Compute elements J+1:N of column J and update the trailing
159 * submatrix.
160 *
161 IF( J.LT.N ) THEN
162 CALL ZDSCAL( N-J, ONE / AJJ, AP( JJ+1 ), 1 )
163 CALL ZHPR( 'Lower', N-J, -ONE, AP( JJ+1 ), 1,
164 $ AP( JJ+N-J+1 ) )
165 JJ = JJ + N - J + 1
166 END IF
167 20 CONTINUE
168 END IF
169 GO TO 40
170 *
171 30 CONTINUE
172 INFO = J
173 *
174 40 CONTINUE
175 RETURN
176 *
177 * End of ZPPTRF
178 *
179 END