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++)
32 Max(double a
, double b
)
34 return(a
> b
? a
: b
);
38 choldecomp(double **a
, int n
, double maxoffl
, double *maxadd
)
40 double minl
, minljj
, minl2
;
43 minl
= pow(macheps(), 0.25) * maxoffl
;
47 for (i
= 0; i
< n
; i
++)
48 maxoffl
= Max(fabs(a
[i
][i
]), maxoffl
);
49 maxoffl
= sqrt(maxoffl
);
50 minl2
= sqrt(macheps()) * maxoffl
;
54 for (j
= 0; j
< n
; j
++) {
55 for (i
= 0; i
< j
; i
++) a
[j
][j
] -= a
[j
][i
] * a
[j
][i
];
59 for (i
= j
+ 1; i
< n
; 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
]);
69 if (minljj
< minl2
) minljj
= minl2
;
70 *maxadd
= Max(*maxadd
, minljj
* minljj
- a
[j
][j
]);
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;