clean up of consistent-dataframe-p generic/methods. They were subtly wrong.
[CommonLispStat.git] / lib / cholesky.c
blob2829d46570382dc705e9cd4b1e9bb8f9f4d026a6
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. */
7 #include "linalg.h"
9 /*
10 choldecomp(a, n)
11 RMatrix a;
12 int n;
14 int i, j, k;
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++)
27 a[i][j] = 0.0;
31 static double
32 Max(double a, double b)
34 return(a > b ? a : b);
37 void
38 choldecomp(double **a, size_t n, double maxoffl, double *maxadd)
40 double minl, minljj, minl2;
41 size_t i, j, k;
43 minl = pow(macheps(), 0.25) * maxoffl;
44 minl2 = 0.0;
46 if (maxoffl == 0.0) {
47 for (i = 0; i < n; i++)
48 maxoffl = Max(fabs(a[i][i]), maxoffl);
49 maxoffl = sqrt(maxoffl);
50 minl2 = sqrt(macheps()) * maxoffl;
53 *maxadd = 0.0;
54 for (j = 0; j < n; j++) {
55 for (i = 0; i < j; i++) a[j][j] -= a[j][i] * a[j][i];
57 minljj = 0.0;
59 for (i = j + 1; i < n; i++) {
60 a[i][j] = a[j][i];
61 for (k = 0; k < j; k++) a[i][j] -= a[i][k] * a[j][k];
62 minljj = Max(fabs(a[i][j]), minljj);
65 minljj = Max(minljj / maxoffl, minl);
67 if (a[j][j] > minljj * minljj) a[j][j] = sqrt(a[j][j]);
68 else {
69 if (minljj < minl2) minljj = minl2;
70 *maxadd = Max(*maxadd, minljj * minljj - a[j][j]);
71 a[j][j] = minljj;
74 for (i = j + 1; i < n; i++) a[i][j] /= a[j][j];
77 for (i = 0; i < n; i++)
78 for (j = i + 1; j < n; j++) a[i][j] = 0.0;