share/linearalgebra/eigens-by-jacobi.lisp

   1 ;;  Copyright 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 (eval-when (:compile-toplevel :load-toplevel :execute)
  11   ($put '$eigensbyjacobi 1 '$version)) ;; Let's have version numbers 1,2,3,...
  12
  13 ;; One sweep zeros each member of the matrix; for a n x n matrix, this requires n(n-1)/2
  14 ;; Jacobi rotations.
  15
  16 ;; For flonum floats, eps is the machine epsilon; for bigfloats, it is 1/2^fpprec.
  17
  18 ;; The variable 'change' tracks the greatest percent change in a diagonal entry in
  19 ;; a sweep.  When the diagonal entry is less than eps, the percent change set to zero.
  20 ;; The iteration stops when 'change' is less than eps (numerical convergence).
  21
  22 ;; The matrix entries are computed with numerically friendly formulae--they
  23 ;; have the form new value <-- old value + correction.  In general, the
  24 ;; correction is 'small.' These formulae are well-known; I used the reference
  25 ;; "Numerical Recipes in Fortran," by Press et.al.
  26
  27 (defun $eigens_by_jacobi (mm &optional (fld-name '$floatfield))
  28   (if (not (member fld-name `($floatfield $bigfloatfield)))
  29       (merror "The field must either be 'floatfield' or 'bigfloatfield'"))
  30
  31   (setq mm (mfuncall '$mat_fullunblocker mm))
  32   ($require_real_symmetric_matrix mm '$first '$eigens_by_jacobi)
  33
  34   (let* ((mat) (g) (h) (sweeps 0) (rotations 0) (eps) (change)
  35          (theta) (mpq) (c) (s)  (tee) (tau) (d) (v) (x) (row)
  36          (n ($first ($matrix_size mm))) (continue (> n 1))
  37          (fld ($require_ring fld-name '$second '$eigens_by_jacobi))
  38          (one (funcall (mring-mult-id fld)))
  39          (zero (funcall (mring-add-id fld))))
  40
  41     (flet
  42         ((fzerop (a) (funcall (mring-fzerop fld) a))
  43          (fabs (a) (funcall (mring-abs fld) a))
  44          (fnegate (a) (funcall (mring-negate fld) a))
  45          (fpsqrt (a) (funcall (mring-psqrt fld) a))
  46          (fadd (a b) (funcall (mring-add fld) a b))
  47          (fsub (a b) (funcall (mring-sub fld) a b))
  48          (fmult (a b) (funcall (mring-mult fld) a b))
  49          (fdiv (a b) (funcall (mring-div fld) a b))
  50          (fgreat (a b) (funcall (mring-great fld) a b))
  51          (fmax (a b) (if (funcall (mring-great fld) a b) a b))
  52          (fconvert (a) (funcall (mring-maxima-to-mring fld) a)))
  53
  54       (setq mat (make-array (list n n) :initial-contents (mapcar #'rest
  55                                                                  (margs (matrix-map #'fconvert mm)))))
  56       (setq v (make-array (list n n) :initial-element zero))
  57       (setq d (make-array n))
  58
  59       (setq eps (if (eq fld-name '$floatfield) flonum-epsilon ($bfloat (div 1 (power 2 fpprec)))))
  60
  61       (decf n)
  62       (loop for i from 0 to n do
  63         (setf (aref v i i) one)
  64         (setf (aref d i) (aref mat i i)))
  65
  66       (while continue
  67         (if (> sweeps 50) (merror "Exceeded maximum allowable number of Jacobi sweeps"))
  68         (incf sweeps)
  69         (loop for p from 0 to n do
  70           (loop for q from (+ p 1) to n do
  71             (setq mpq (aref mat p q))
  72             (cond ((not (fzerop mpq))
  73                    (incf rotations)
  74                    (setq theta (fdiv (fsub (aref mat q q) (aref mat p p))(fmult 2 mpq)))
  75                    (setq tee (fdiv one (fadd (fabs theta) (fpsqrt (fadd one (fmult theta theta))))))
  76                    (if (fgreat 0 theta) (setq tee (fnegate tee)))
  77                    (setq c (fdiv one (fpsqrt (fadd one (fmult tee tee)))))
  78                    (setq s (fmult tee c))
  79                    (setq tau (fdiv s (fadd one c)))
  80                    (setf (aref mat p q) zero)
  81
  82                    (loop for k from 0 to (- p 1) do
  83                      (setq g (aref mat k p))
  84                      (setq h (aref mat k q))
  85                      (setf (aref mat k p) (fsub g (fmult s (fadd h (fmult g tau)))))
  86                      (setf (aref mat k q) (fadd h (fmult s (fsub g (fmult h tau))))))
  87
  88                    (loop for k from (+ p 1) to (- q 1) do
  89                      (setq g (aref mat p k))
  90                      (setq h (aref mat k q))
  91                      (setf (aref mat p k) (fsub g (fmult s (fadd h (fmult g tau)))))
  92                      (setf (aref mat k q) (fadd h (fmult s (fsub g (fmult h tau))))))
  93
  94                    (loop for k from (+ q 1) to n do
  95                      (setq g (aref mat p k))
  96                      (setq h (aref mat q k))
  97                      (setf (aref mat p k) (fsub g (fmult s (fadd h (fmult g tau)))))
  98                      (setf (aref mat q k) (fadd h (fmult s (fsub g (fmult h tau))))))
  99
 100                    (setf (aref mat p p) (fsub (aref mat p p) (fmult tee mpq)))
 101                    (setf (aref mat q q) (fadd (aref mat q q) (fmult tee mpq)))
 102                    (loop for k from 0 to n do
 103                      (setq g (aref v k p))
 104                      (setq h (aref v k q))
 105                      (setf (aref v k p) (fsub g (fmult s (fadd h (fmult g tau)))))
 106                      (setf (aref v k q)(fadd h (fmult s (fsub g (fmult h tau))))))))))
 107
 108         (setq change zero)
 109         (loop for i from 0 to n do
 110           (setq x (aref mat i i))
 111           (setq change (fmax change (if (fgreat (fabs x) eps) (fabs (fdiv (fsub (aref d i) x) x)) zero)))
 112           (setf (aref d i) x))
 113
 114         (inform '$debug '$linearalgebra "The largest percent change was ~:M~%" change)
 115         (setq continue (fgreat change eps)))
 116
 117       (inform '$verbose '$linearalgebra "number of sweeps: ~:M~%" sweeps)
 118       (inform '$verbose '$linearalgebra "number of rotations: ~:M~%" rotations)
 119
 120       (setq mm nil)
 121       (loop for i from 0 to n do
 122         (setq row nil)
 123         (loop for j from 0 to n do
 124           (push (aref v i j) row))
 125         (setq row (reverse row))
 126         (push '(mlist) row)
 127         (push row mm))
 128       (setq mm (reverse mm))
 129       (push '($matrix) mm)
 130       (setq d `((mlist) ,@(coerce d 'list)))
 131       `((mlist) ,d ,mm))))
 132
 133
 134
 135
 136
 137
 138
 139
 140
 141
 142
 143
 144
 145