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
  32 Max(double a, double b)
  33 {
  34   return(a > b ? a : b);
  35 }
  36
  37 void
  38 choldecomp(double **a, size_t n, double maxoffl, double *maxadd)
  39 {
  40   double minl, minljj, minl2;
  41   size_t i, j, k;
  42
  43   minl = pow(macheps(), 0.25) * maxoffl;
  44   minl2 = 0.0;
  45
  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;
  51   }
  52
  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];
  56
  57     minljj = 0.0;
  58
  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);
  63     }
  64
  65     minljj = Max(minljj / maxoffl, minl);
  66
  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;
  72     }
  73
  74     for (i = j + 1; i < n; i++) a[i][j] /= a[j][j];
  75   }
  76
  77   for (i = 0; i < n; i++)
  78     for (j = i + 1; j < n; j++) a[i][j] = 0.0;
  79 }