1 /* ./src_f77/cpptrf.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
8 /* Table of constant values */
9
10 static integer c__1 = 1;
11 static real c_b16 = -1.f;
12
cpptrf_(char * uplo,integer * n,complex * ap,integer * info,ftnlen uplo_len)13 /* Subroutine */ int cpptrf_(char *uplo, integer *n, complex *ap, integer *
14 info, ftnlen uplo_len)
15 {
16 /* System generated locals */
17 integer i__1, i__2, i__3;
18 real r__1;
19 complex q__1, q__2;
20
21 /* Builtin functions */
22 double sqrt(doublereal);
23
24 /* Local variables */
25 static integer j, jc, jj;
26 static real ajj;
27 extern /* Subroutine */ int chpr_(char *, integer *, real *, complex *,
28 integer *, complex *, ftnlen);
29 extern /* Complex */ VOID cdotc_(complex *, integer *, complex *, integer
30 *, complex *, integer *);
31 extern logical lsame_(char *, char *, ftnlen, ftnlen);
32 static logical upper;
33 extern /* Subroutine */ int ctpsv_(char *, char *, char *, integer *,
34 complex *, complex *, integer *, ftnlen, ftnlen, ftnlen), csscal_(
35 integer *, real *, complex *, integer *), xerbla_(char *, integer
36 *, ftnlen);
37
38
39 /* -- LAPACK routine (version 3.0) -- */
40 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
41 /* Courant Institute, Argonne National Lab, and Rice University */
42 /* September 30, 1994 */
43
44 /* .. Scalar Arguments .. */
45 /* .. */
46 /* .. Array Arguments .. */
47 /* .. */
48
49 /* Purpose */
50 /* ======= */
51
52 /* CPPTRF computes the Cholesky factorization of a complex Hermitian */
53 /* positive definite matrix A stored in packed format. */
54
55 /* The factorization has the form */
56 /* A = U**H * U, if UPLO = 'U', or */
57 /* A = L * L**H, if UPLO = 'L', */
58 /* where U is an upper triangular matrix and L is lower triangular. */
59
60 /* Arguments */
61 /* ========= */
62
63 /* UPLO (input) CHARACTER*1 */
64 /* = 'U': Upper triangle of A is stored; */
65 /* = 'L': Lower triangle of A is stored. */
66
67 /* N (input) INTEGER */
68 /* The order of the matrix A. N >= 0. */
69
70 /* AP (input/output) COMPLEX array, dimension (N*(N+1)/2) */
71 /* On entry, the upper or lower triangle of the Hermitian matrix */
72 /* A, packed columnwise in a linear array. The j-th column of A */
73 /* is stored in the array AP as follows: */
74 /* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; */
75 /* if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. */
76 /* See below for further details. */
77
78 /* On exit, if INFO = 0, the triangular factor U or L from the */
79 /* Cholesky factorization A = U**H*U or A = L*L**H, in the same */
80 /* storage format as A. */
81
82 /* INFO (output) INTEGER */
83 /* = 0: successful exit */
84 /* < 0: if INFO = -i, the i-th argument had an illegal value */
85 /* > 0: if INFO = i, the leading minor of order i is not */
86 /* positive definite, and the factorization could not be */
87 /* completed. */
88
89 /* Further Details */
90 /* =============== */
91
92 /* The packed storage scheme is illustrated by the following example */
93 /* when N = 4, UPLO = 'U': */
94
95 /* Two-dimensional storage of the Hermitian matrix A: */
96
97 /* a11 a12 a13 a14 */
98 /* a22 a23 a24 */
99 /* a33 a34 (aij = conjg(aji)) */
100 /* a44 */
101
102 /* Packed storage of the upper triangle of A: */
103
104 /* AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ] */
105
106 /* ===================================================================== */
107
108 /* .. Parameters .. */
109 /* .. */
110 /* .. Local Scalars .. */
111 /* .. */
112 /* .. External Functions .. */
113 /* .. */
114 /* .. External Subroutines .. */
115 /* .. */
116 /* .. Intrinsic Functions .. */
117 /* .. */
118 /* .. Executable Statements .. */
119
120 /* Test the input parameters. */
121
122 /* Parameter adjustments */
123 --ap;
124
125 /* Function Body */
126 *info = 0;
127 upper = lsame_(uplo, "U", (ftnlen)1, (ftnlen)1);
128 if (! upper && ! lsame_(uplo, "L", (ftnlen)1, (ftnlen)1)) {
129 *info = -1;
130 } else if (*n < 0) {
131 *info = -2;
132 }
133 if (*info != 0) {
134 i__1 = -(*info);
135 xerbla_("CPPTRF", &i__1, (ftnlen)6);
136 return 0;
137 }
138
139 /* Quick return if possible */
140
141 if (*n == 0) {
142 return 0;
143 }
144
145 if (upper) {
146
147 /* Compute the Cholesky factorization A = U'*U. */
148
149 jj = 0;
150 i__1 = *n;
151 for (j = 1; j <= i__1; ++j) {
152 jc = jj + 1;
153 jj += j;
154
155 /* Compute elements 1:J-1 of column J. */
156
157 if (j > 1) {
158 i__2 = j - 1;
159 ctpsv_("Upper", "Conjugate transpose", "Non-unit", &i__2, &ap[
160 1], &ap[jc], &c__1, (ftnlen)5, (ftnlen)19, (ftnlen)8);
161 }
162
163 /* Compute U(J,J) and test for non-positive-definiteness. */
164
165 i__2 = jj;
166 r__1 = ap[i__2].r;
167 i__3 = j - 1;
168 cdotc_(&q__2, &i__3, &ap[jc], &c__1, &ap[jc], &c__1);
169 q__1.r = r__1 - q__2.r, q__1.i = -q__2.i;
170 ajj = q__1.r;
171 if (ajj <= 0.f) {
172 i__2 = jj;
173 ap[i__2].r = ajj, ap[i__2].i = 0.f;
174 goto L30;
175 }
176 i__2 = jj;
177 r__1 = sqrt(ajj);
178 ap[i__2].r = r__1, ap[i__2].i = 0.f;
179 /* L10: */
180 }
181 } else {
182
183 /* Compute the Cholesky factorization A = L*L'. */
184
185 jj = 1;
186 i__1 = *n;
187 for (j = 1; j <= i__1; ++j) {
188
189 /* Compute L(J,J) and test for non-positive-definiteness. */
190
191 i__2 = jj;
192 ajj = ap[i__2].r;
193 if (ajj <= 0.f) {
194 i__2 = jj;
195 ap[i__2].r = ajj, ap[i__2].i = 0.f;
196 goto L30;
197 }
198 ajj = sqrt(ajj);
199 i__2 = jj;
200 ap[i__2].r = ajj, ap[i__2].i = 0.f;
201
202 /* Compute elements J+1:N of column J and update the trailing */
203 /* submatrix. */
204
205 if (j < *n) {
206 i__2 = *n - j;
207 r__1 = 1.f / ajj;
208 csscal_(&i__2, &r__1, &ap[jj + 1], &c__1);
209 i__2 = *n - j;
210 chpr_("Lower", &i__2, &c_b16, &ap[jj + 1], &c__1, &ap[jj + *n
211 - j + 1], (ftnlen)5);
212 jj = jj + *n - j + 1;
213 }
214 /* L20: */
215 }
216 }
217 goto L40;
218
219 L30:
220 *info = j;
221
222 L40:
223 return 0;
224
225 /* End of CPPTRF */
226
227 } /* cpptrf_ */
228
229