1 /* ./src_f77/spotf2.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_b10 = -1.f;
12 static real c_b12 = 1.f;
13
spotf2_(char * uplo,integer * n,real * a,integer * lda,integer * info,ftnlen uplo_len)14 /* Subroutine */ int spotf2_(char *uplo, integer *n, real *a, integer *lda,
15 integer *info, ftnlen uplo_len)
16 {
17 /* System generated locals */
18 integer a_dim1, a_offset, i__1, i__2, i__3;
19 real r__1;
20
21 /* Builtin functions */
22 double sqrt(doublereal);
23
24 /* Local variables */
25 static integer j;
26 static real ajj;
27 extern doublereal sdot_(integer *, real *, integer *, real *, integer *);
28 extern logical lsame_(char *, char *, ftnlen, ftnlen);
29 extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *),
30 sgemv_(char *, integer *, integer *, real *, real *, integer *,
31 real *, integer *, real *, real *, integer *, ftnlen);
32 static logical upper;
33 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
34
35
36 /* -- LAPACK routine (version 3.0) -- */
37 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
38 /* Courant Institute, Argonne National Lab, and Rice University */
39 /* February 29, 1992 */
40
41 /* .. Scalar Arguments .. */
42 /* .. */
43 /* .. Array Arguments .. */
44 /* .. */
45
46 /* Purpose */
47 /* ======= */
48
49 /* SPOTF2 computes the Cholesky factorization of a real symmetric */
50 /* positive definite matrix A. */
51
52 /* The factorization has the form */
53 /* A = U' * U , if UPLO = 'U', or */
54 /* A = L * L', if UPLO = 'L', */
55 /* where U is an upper triangular matrix and L is lower triangular. */
56
57 /* This is the unblocked version of the algorithm, calling Level 2 BLAS. */
58
59 /* Arguments */
60 /* ========= */
61
62 /* UPLO (input) CHARACTER*1 */
63 /* Specifies whether the upper or lower triangular part of the */
64 /* symmetric matrix A is stored. */
65 /* = 'U': Upper triangular */
66 /* = 'L': Lower triangular */
67
68 /* N (input) INTEGER */
69 /* The order of the matrix A. N >= 0. */
70
71 /* A (input/output) REAL array, dimension (LDA,N) */
72 /* On entry, the symmetric matrix A. If UPLO = 'U', the leading */
73 /* n by n upper triangular part of A contains the upper */
74 /* triangular part of the matrix A, and the strictly lower */
75 /* triangular part of A is not referenced. If UPLO = 'L', the */
76 /* leading n by n lower triangular part of A contains the lower */
77 /* triangular part of the matrix A, and the strictly upper */
78 /* triangular part of A is not referenced. */
79
80 /* On exit, if INFO = 0, the factor U or L from the Cholesky */
81 /* factorization A = U'*U or A = L*L'. */
82
83 /* LDA (input) INTEGER */
84 /* The leading dimension of the array A. LDA >= max(1,N). */
85
86 /* INFO (output) INTEGER */
87 /* = 0: successful exit */
88 /* < 0: if INFO = -k, the k-th argument had an illegal value */
89 /* > 0: if INFO = k, the leading minor of order k is not */
90 /* positive definite, and the factorization could not be */
91 /* completed. */
92
93 /* ===================================================================== */
94
95 /* .. Parameters .. */
96 /* .. */
97 /* .. Local Scalars .. */
98 /* .. */
99 /* .. External Functions .. */
100 /* .. */
101 /* .. External Subroutines .. */
102 /* .. */
103 /* .. Intrinsic Functions .. */
104 /* .. */
105 /* .. Executable Statements .. */
106
107 /* Test the input parameters. */
108
109 /* Parameter adjustments */
110 a_dim1 = *lda;
111 a_offset = 1 + a_dim1;
112 a -= a_offset;
113
114 /* Function Body */
115 *info = 0;
116 upper = lsame_(uplo, "U", (ftnlen)1, (ftnlen)1);
117 if (! upper && ! lsame_(uplo, "L", (ftnlen)1, (ftnlen)1)) {
118 *info = -1;
119 } else if (*n < 0) {
120 *info = -2;
121 } else if (*lda < max(1,*n)) {
122 *info = -4;
123 }
124 if (*info != 0) {
125 i__1 = -(*info);
126 xerbla_("SPOTF2", &i__1, (ftnlen)6);
127 return 0;
128 }
129
130 /* Quick return if possible */
131
132 if (*n == 0) {
133 return 0;
134 }
135
136 if (upper) {
137
138 /* Compute the Cholesky factorization A = U'*U. */
139
140 i__1 = *n;
141 for (j = 1; j <= i__1; ++j) {
142
143 /* Compute U(J,J) and test for non-positive-definiteness. */
144
145 i__2 = j - 1;
146 ajj = a[j + j * a_dim1] - sdot_(&i__2, &a[j * a_dim1 + 1], &c__1,
147 &a[j * a_dim1 + 1], &c__1);
148 if (ajj <= 0.f) {
149 a[j + j * a_dim1] = ajj;
150 goto L30;
151 }
152 ajj = sqrt(ajj);
153 a[j + j * a_dim1] = ajj;
154
155 /* Compute elements J+1:N of row J. */
156
157 if (j < *n) {
158 i__2 = j - 1;
159 i__3 = *n - j;
160 sgemv_("Transpose", &i__2, &i__3, &c_b10, &a[(j + 1) * a_dim1
161 + 1], lda, &a[j * a_dim1 + 1], &c__1, &c_b12, &a[j + (
162 j + 1) * a_dim1], lda, (ftnlen)9);
163 i__2 = *n - j;
164 r__1 = 1.f / ajj;
165 sscal_(&i__2, &r__1, &a[j + (j + 1) * a_dim1], lda);
166 }
167 /* L10: */
168 }
169 } else {
170
171 /* Compute the Cholesky factorization A = L*L'. */
172
173 i__1 = *n;
174 for (j = 1; j <= i__1; ++j) {
175
176 /* Compute L(J,J) and test for non-positive-definiteness. */
177
178 i__2 = j - 1;
179 ajj = a[j + j * a_dim1] - sdot_(&i__2, &a[j + a_dim1], lda, &a[j
180 + a_dim1], lda);
181 if (ajj <= 0.f) {
182 a[j + j * a_dim1] = ajj;
183 goto L30;
184 }
185 ajj = sqrt(ajj);
186 a[j + j * a_dim1] = ajj;
187
188 /* Compute elements J+1:N of column J. */
189
190 if (j < *n) {
191 i__2 = *n - j;
192 i__3 = j - 1;
193 sgemv_("No transpose", &i__2, &i__3, &c_b10, &a[j + 1 +
194 a_dim1], lda, &a[j + a_dim1], lda, &c_b12, &a[j + 1 +
195 j * a_dim1], &c__1, (ftnlen)12);
196 i__2 = *n - j;
197 r__1 = 1.f / ajj;
198 sscal_(&i__2, &r__1, &a[j + 1 + j * a_dim1], &c__1);
199 }
200 /* L20: */
201 }
202 }
203 goto L40;
204
205 L30:
206 *info = j;
207
208 L40:
209 return 0;
210
211 /* End of SPOTF2 */
212
213 } /* spotf2_ */
214
215