1 /* choldecomp - Cholesky decomposition routines. */
2 /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney */
3 /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz */
4 /* You may give out copies of this software; for conditions see the */
5 /* file COPYING included with this distribution. */
16 for (i = 0; i < n; i++) {
17 for (k = 0; k < i; k++)
18 a[i][i] -= a[k][i] * a[k][i];
19 a[i][i] = (a[i][i] > 0) ? sqrt(a[i][i]) : 0.0;
20 for (j = i + 1; j < n; j++) {
21 for (k = 0; k < i; k++) a[i][j] -= a[k][i] * a[k][j];
22 a[i][j] = (a[i][i] > 0.0) ? a[i][j] / a[i][i] : 0.0;
25 for (i = 0; i < n; i++)
26 for (j = 0; j < i; j++)
31 static double Max(a
, b
)
34 return(a
> b
? a
: b
);
37 choldecomp(a
, n
, maxoffl
, maxadd
)
40 double maxoffl
, *maxadd
;
42 double minl
, minljj
, minl2
;
45 minl
= pow(macheps(), 0.25) * maxoffl
;
49 for (i
= 0; i
< n
; i
++)
50 maxoffl
= Max(fabs(a
[i
][i
]), maxoffl
);
51 maxoffl
= sqrt(maxoffl
);
52 minl2
= sqrt(macheps()) * maxoffl
;
56 for (j
= 0; j
< n
; j
++) {
57 for (i
= 0; i
< j
; i
++) a
[j
][j
] -= a
[j
][i
] * a
[j
][i
];
61 for (i
= j
+ 1; i
< n
; i
++) {
63 for (k
= 0; k
< j
; k
++) a
[i
][j
] -= a
[i
][k
] * a
[j
][k
];
64 minljj
= Max(fabs(a
[i
][j
]), minljj
);
67 minljj
= Max(minljj
/ maxoffl
, minl
);
69 if (a
[j
][j
] > minljj
* minljj
) a
[j
][j
] = sqrt(a
[j
][j
]);
71 if (minljj
< minl2
) minljj
= minl2
;
72 *maxadd
= Max(*maxadd
, minljj
* minljj
- a
[j
][j
]);
76 for (i
= j
+ 1; i
< n
; i
++) a
[i
][j
] /= a
[j
][j
];
79 for (i
= 0; i
< n
; i
++)
80 for (j
= i
+ 1; j
< n
; j
++) a
[i
][j
] = 0.0;