share/linearalgebra/linalgcholesky.lisp

   1 ;;  Copyright 2005, 2006 by Barton Willis
   2
   3 ;;  This is free software; you can redistribute it and/or
   4 ;;  modify it under the terms of the GNU General Public License,
   5 ;;  http://www.gnu.org/copyleft/gpl.html.
   6
   7 ;; This software has NO WARRANTY, not even the implied warranty of
   8 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
   9
  10 ($put '$cholesky 1 '$version)
  11         (defun $cholesky (m &optional (fld-name '$generalring))
  12   ($require_nonempty_matrix m '$first '$cholesky)
  13   ($require_selfadjoint_matrix m '$first '$cholesky)
  14
  15   (let* ((n ($first ($matrix_size m))) (lii) (lii-inv) (l) (acc) (mat) (row)
  16          (fld ($require_ring fld-name '$second '$cholesky)))
  17          ;;(zero (funcall (mring-add-id fld))))
  18
  19     (flet
  20         ((fzerop (a) (funcall (mring-fzerop fld) a))
  21          (fpsqrt (a) (funcall (mring-psqrt fld) a))
  22          (fsub (a b) (funcall (mring-sub fld) a b))
  23          (fmult (a b) (funcall (mring-mult fld) a b))
  24          (fadjoint (a) (funcall (mring-adjoint fld) a))
  25          (freciprocal (a) (funcall (mring-reciprocal fld) a))
  26          (frevert (a) (funcall (mring-mring-to-maxima fld) a))
  27          (fconvert (a) (funcall (mring-maxima-to-mring fld) a)))
  28
  29       (setq l ($zerofor m))
  30       (setq l (make-array (list n n) :initial-contents (mapcar #'rest (margs l))))
  31       (setq mat (make-array (list n n) :initial-contents (mapcar #'rest (margs m))))
  32
  33       (decf n)
  34
  35       ;; Convert each entry of mat to a ring member.
  36
  37       (loop for i from 0 to n do
  38         (loop for j from 0 to n do
  39           (setf (aref mat i j) (fconvert (aref mat i j)))))
  40
  41
  42       (loop for i from 0 to n do
  43         (setq acc (aref mat i i))
  44         (loop for k from 0 to (- i 1) do
  45           (setq acc (fsub acc (fmult (aref l i k) (fadjoint (aref l i k))))))
  46
  47         (setq lii (if ($matrixp acc) ($cholesky acc fld-name) (fpsqrt acc)))
  48         (if (null lii) (merror "Unable to find the Cholesky factorization"))
  49         (setf (aref l i i) lii)
  50         (cond ((<= (+ i 1) n)
  51                (if (fzerop lii) (merror "Unable to find the Cholesky factorization"))
  52                (setq lii-inv (fadjoint (freciprocal lii)))))
  53         (loop for j from (+ i 1) to n do
  54           (setq acc (aref mat j i))
  55           (loop for k from 0 to (- i 1) do
  56             (setq acc (fsub acc (fmult (aref l j k) (fadjoint (aref l i k))))))
  57           (setf (aref l j i) (fmult acc lii-inv))))
  58
  59       (setq mat nil)
  60       (loop for i from 0 to n do
  61         (setq row nil)
  62         (loop for j from 0 to n do
  63           (push (full-matrix-map  #'frevert (aref l i j)) row))
  64         (setq row (reverse row))
  65         (push '(mlist) row)
  66         (push row mat))
  67       (setq mat (reverse mat))
  68       (push '($matrix) mat))))
  69
  70