lib/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 "linalg.h"
   8
   9 /*
  10 choldecomp(a, n)
  11         RMatrix a;
  12         int n;
  13 {
  14   int i, j, k;
  15
  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;
  23     }
  24   }
  25   for (i = 0; i < n; i++)
  26     for (j = 0; j < i; j++)
  27       a[i][j] = 0.0;
  28 }
  29 */
  30
  31 static double Max(a, b)
  32         double a, b;
  33 {
  34   return(a > b ? a : b);
  35 }
  36
  37 choldecomp(a, n, maxoffl, maxadd)
  38         RMatrix a;
  39         int n;
  40         double maxoffl, *maxadd;
  41 {
  42   double minl, minljj, minl2;
  43   int i, j, k;
  44
  45   minl = pow(macheps(), 0.25) * maxoffl;
  46   minl2 = 0.0;
  47
  48   if (maxoffl == 0.0) {
  49     for (i = 0; i < n; i++)
  50       maxoffl = Max(fabs(a[i][i]), maxoffl);
  51     maxoffl = sqrt(maxoffl);
  52     minl2 = sqrt(macheps()) * maxoffl;
  53   }
  54
  55   *maxadd = 0.0;
  56   for (j = 0; j < n; j++) {
  57     for (i = 0; i < j; i++) a[j][j] -= a[j][i] * a[j][i];
  58
  59     minljj = 0.0;
  60
  61     for (i = j + 1; i < n; i++) {
  62       a[i][j] = a[j][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);
  65     }
  66
  67     minljj = Max(minljj / maxoffl, minl);
  68
  69     if (a[j][j] > minljj * minljj) a[j][j] = sqrt(a[j][j]);
  70     else {
  71       if (minljj < minl2) minljj = minl2;
  72       *maxadd = Max(*maxadd, minljj * minljj - a[j][j]);
  73       a[j][j] = minljj;
  74     }
  75
  76     for (i = j + 1; i < n; i++) a[i][j] /= a[j][j];
  77   }
  78
  79   for (i = 0; i < n; i++)
  80     for (j = i + 1; j < n; j++) a[i][j] = 0.0;
  81 }