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
146