1 /* dl7srt.f -- translated by f2c (version 19970211).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include "f2c.h"
7
dl7srt_(n1,n,l,a,irc)8 /* Subroutine */ int dl7srt_(n1, n, l, a, irc)
9 integer *n1, *n;
10 doublereal *l, *a;
11 integer *irc;
12 {
13 /* System generated locals */
14 integer i__1, i__2, i__3;
15
16 /* Builtin functions */
17 double sqrt();
18
19 /* Local variables */
20 static integer i__, j, k;
21 static doublereal t;
22 static integer i0, j0, ij, ik, jk;
23 static doublereal td;
24 static integer im1, jm1;
25
26
27 /* *** COMPUTE ROWS N1 THROUGH N OF THE CHOLESKY FACTOR L OF */
28 /* *** A = L*(L**T), WHERE L AND THE LOWER TRIANGLE OF A ARE BOTH
29 */
30 /* *** STORED COMPACTLY BY ROWS (AND MAY OCCUPY THE SAME STORAGE). */
31 /* *** IRC = 0 MEANS ALL WENT WELL. IRC = J MEANS THE LEADING */
32 /* *** PRINCIPAL J X J SUBMATRIX OF A IS NOT POSITIVE DEFINITE -- */
33 /* *** AND L(J*(J+1)/2) CONTAINS THE (NONPOS.) REDUCED J-TH DIAGONAL.
34 */
35
36 /* *** PARAMETERS *** */
37
38 /* DIMENSION L(N*(N+1)/2), A(N*(N+1)/2) */
39
40 /* *** LOCAL VARIABLES *** */
41
42
43 /* *** INTRINSIC FUNCTIONS *** */
44 /* /+ */
45 /* / */
46 /* /6 */
47 /* DATA ZERO/0.D+0/ */
48 /* /7 */
49 /* / */
50
51 /* *** BODY *** */
52
53 /* Parameter adjustments */
54 --a;
55 --l;
56
57 /* Function Body */
58 i0 = *n1 * (*n1 - 1) / 2;
59 i__1 = *n;
60 for (i__ = *n1; i__ <= i__1; ++i__) {
61 td = 0.;
62 if (i__ == 1) {
63 goto L40;
64 }
65 j0 = 0;
66 im1 = i__ - 1;
67 i__2 = im1;
68 for (j = 1; j <= i__2; ++j) {
69 t = 0.;
70 if (j == 1) {
71 goto L20;
72 }
73 jm1 = j - 1;
74 i__3 = jm1;
75 for (k = 1; k <= i__3; ++k) {
76 ik = i0 + k;
77 jk = j0 + k;
78 t += l[ik] * l[jk];
79 /* L10: */
80 }
81 L20:
82 ij = i0 + j;
83 j0 += j;
84 t = (a[ij] - t) / l[j0];
85 l[ij] = t;
86 td += t * t;
87 /* L30: */
88 }
89 L40:
90 i0 += i__;
91 t = a[i0] - td;
92 if (t <= 0.) {
93 goto L60;
94 }
95 l[i0] = sqrt(t);
96 /* L50: */
97 }
98
99 *irc = 0;
100 goto L999;
101
102 L60:
103 l[i0] = t;
104 *irc = i__;
105
106 L999:
107 return 0;
108
109 /* *** LAST CARD OF DL7SRT *** */
110 } /* dl7srt_ */
111
112