src/c/cholesky.c

   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.                       */
   6
   7 #include "xlisp.h"
   8 #include "xlstat.h"
   9 #include "linalg.h"
  10
  11 /* forward declarations */
  12 LOCAL double Max P2H(double, double);
  13
  14 LOCAL double Max P2C(double, a, double, b)
  15 {
  16   return(a > b ? a : b);
  17 }
  18
  19 #define Dot(n, x, y) blas_ddot(n, x, 1, y, 1)
  20
  21 VOID choldecomp P4C(double *, a, int, n, double, maxoffl, double *, maxadd)
  22 {
  23   double minl, minljj, minl2;
  24   int i, j;
  25   int in, jn;
  26
  27   minl = pow(macheps(), 0.25) * maxoffl;
  28   minl2 = 0.0;
  29
  30   if (maxoffl == 0.0) {
  31     for (i = 0; i < n; i++) {
  32       int ii = i * n + i;
  33       maxoffl = Max(fabs(a[ii]), maxoffl);
  34     }
  35     maxoffl = sqrt(maxoffl);
  36     minl2 = sqrt(macheps()) * maxoffl;
  37   }
  38
  39   *maxadd = 0.0;
  40   for (j = 0, jn = 0; j < n; j++, jn += n) {
  41     int jj = jn + j;
  42     a[jj] -= Dot(j, a + jn, a + jn);
  43
  44     minljj = 0.0;
  45
  46     for (i = j + 1, in = (j + 1) * n; i < n; i++, in += n) {
  47       int ij = in + j;
  48       int ji = jn + i;
  49       a[ij] = a[ji];
  50       a[ij] -= Dot(j, a + in, a + jn);
  51       minljj = Max(fabs(a[ij]), minljj);
  52     }
  53
  54     minljj = Max(minljj / maxoffl, minl);
  55
  56     if (a[jj] > minljj * minljj) a[jj] = sqrt(a[jj]);
  57     else {
  58       if (minljj < minl2) minljj = minl2;
  59       *maxadd = Max(*maxadd, minljj * minljj - a[jj]);
  60       a[jj] = minljj;
  61     }
  62
  63     /**** use BLAS dscal */
  64     for (i = j + 1, in = (j + 1) * n; i < n; i++, in += n) {
  65       int ij = in + j;
  66       a[ij] /= a[jj];
  67     }
  68   }
  69
  70   for (i = 0, in = 0; i < n; i++, in += n) {
  71     for (j = i + 1; j < n; j++) {
  72       int ij = in + j;
  73       a[ij] = 0.0;
  74     }
  75   }
  76 }