archive/share/trash/cholesky.mac

   1 /* Compute Cholesky decomposition of A,
   2  * a lower-triangular matrix L such that L . transpose(L) = A
   3  *
   4  * Examples:
   5
   6    A : matrix ([a, b, c], [d, e, f], [g, h, i]);
   7    A2 : transpose (A) . A;
   8    B : cholesky (A2);
   9    B . transpose (B) - A2;
  10
  11    A : matrix ([2, 3, 4], [-2, 2,- 3], [11, -2, 3]);
  12    A2 : transpose (A) . A;
  13    B : cholesky (A2);
  14    B . transpose (B) - A2;
  15
  16  *
  17  * copyright Robert Dodier, 2005/11/01
  18  * Released under the terms of the GNU Public License
  19  *
  20  */
  21 cholesky (A):= block
  22
  23    ([n : length (A),
  24     L : copymatrix (A),
  25     p : makelist (0, i, 1, length (A))],
  26
  27     for i thru n do
  28         for j : i thru n do
  29            (x : L[i, j],
  30             x : x - sum (L[j, k] * L[i, k], k, 1, i - 1),
  31             if i = j then
  32                 p[i] : 1 / sqrt(x)
  33             else
  34                 L[j, i] : x * p[i]),
  35
  36     for i thru n do
  37         L[i, i] : 1 / p[i],
  38
  39     for i thru n do
  40         for j : i + 1 thru n do
  41             L[i, j] : 0,
  42
  43     L);