Skip to content

Commit b264705

Browse files
stylewarningstylewarning
committed
compute eigenvalues (not vectors) in pure Lisp
1 parent f03e967 commit b264705

File tree

1 file changed

+178
-0
lines changed
  • src/high-level/matrix-functions/eig

1 file changed

+178
-0
lines changed

‎src/high-level/matrix-functions/eig/eig.lisp

Lines changed: 178 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,178 @@
1+
(in-package #:magicl)
2+
3+
(defvar *junk-tol* 1d-8)
4+
5+
(defgeneric to-row-major-lisp-array (m)
6+
(:method ((m matrix/double-float))
7+
(destructuring-bind (rows cols) (shape m)
8+
(let ((x (make-array (* rows cols) :element-type 'double-float :initial-element 0.0d0)))
9+
(dotimes (r rows x)
10+
(dotimes (c cols)
11+
(setf (aref x (+ c (* r cols))) (tref m r c))))))))
12+
13+
(defun cplx (a b)
14+
(if (zerop b) a (complex a b)))
15+
16+
(defun internal-eig (m)
17+
(let* ((n (isqrt (length m)))
18+
(eigs-real (make-array n :element-type 'double-float :initial-element 0.0d0))
19+
(eigs-imag (make-array n :element-type 'double-float :initial-element 0.0d0))
20+
(left-vecs (make-array 0 :element-type 'double-float :initial-element 0.0d0))
21+
(right-vecs (make-array (* 2 n n) :element-type 'double-float :initial-element 0.0d0))
22+
(lwork (* 4 n n))
23+
(work (make-array lwork :element-type 'double-float :initial-element 0.0d0)))
24+
(lapack::dgeev
25+
"N" ; left eigenvectors
26+
"V" ; right eigenvectors
27+
n ; order
28+
m ; matrix (flattened)
29+
n ; leading dimension
30+
eigs-real ; eigenvalues (real part)
31+
eigs-imag ; eigenvalues (imag part)
32+
left-vecs ; left eigenvectors (not computed)
33+
(* 2 n)
34+
right-vecs ; right eigenvectors
35+
(* 2 n)
36+
work
37+
lwork
38+
0 ; info
39+
)
40+
(values eigs-real eigs-imag right-vecs)))
41+
42+
(defmethod eig-lisp ((m matrix/double-float))
43+
(assert (square-matrix-p m))
44+
(let* ((a (to-row-major-lisp-array m))
45+
(shape (shape m))
46+
(n (first shape)))
47+
(multiple-value-bind (val-re val-im vecs) (internal-eig a)
48+
(let ((eigenvalues (cl:map 'list #'cplx val-re val-im))
49+
(eigenvectors (zeros shape :type '(complex double-float))))
50+
(loop :with j := 0
51+
:with eigenvalues-left := eigenvalues
52+
:until (null eigenvalues-left)
53+
:for e := (pop eigenvalues-left)
54+
:do (etypecase e
55+
(real
56+
(dotimes (i n)
57+
(setf (tref eigenvectors i j)
58+
(complex (aref vecs (+ j (* i (* 2 n))))
59+
0.0d0)))
60+
(incf j 1))
61+
(complex
62+
(let ((next (pop eigenvalues-left)))
63+
(assert (cl:= next (conjugate e))
64+
()
65+
()
66+
"Expected eigenvalues to come in conjugate pairs. Got ~A then ~A, which don't appear to be conjugates." e next))
67+
(dotimes (i n)
68+
(let ((re (aref vecs (+ j (* i (* 2 n)))))
69+
(im (aref vecs (+ (+ j 1) (* i (* 2 n))))))
70+
(setf (tref eigenvectors i j)
71+
(complex re im)
72+
(tref eigenvectors i (+ j 1))
73+
(complex re (- im)))))
74+
(incf j 2))))
75+
(values eigenvalues
76+
eigenvectors)))))
77+
78+
(defun embed-complex (m)
79+
(assert (square-matrix-p m))
80+
(let* ((n (nrows m))
81+
(embedding (zeros (list (* 2 n) (* 2 n)) :type 'double-float)))
82+
;; map a+bi -> [a -b; b a] for all elements of M
83+
(uiop:nest
84+
(dotimes (z-row n))
85+
(let ((r-row (* 2 z-row))))
86+
(dotimes (z-col n))
87+
(let ((r-col (* 2 z-col))))
88+
(let* ((z (tref m z-row z-col))
89+
(re-z (realpart z))
90+
(im-z (imagpart z))))
91+
(progn
92+
(setf (tref embedding r-row r-col)
93+
re-z
94+
(tref embedding r-row (1+ r-col))
95+
(- im-z)
96+
(tref embedding (1+ r-row) r-col)
97+
im-z
98+
(tref embedding (1+ r-row) (1+ r-col))
99+
re-z)))
100+
embedding))
101+
102+
(defmethod eig-lisp ((m matrix/complex-double-float))
103+
(assert (square-matrix-p m))
104+
(multiple-value-bind (evals evecs)
105+
(eig-lisp (embed-complex m))
106+
(let* ((tr (trace m))
107+
(tr-real (realpart tr))
108+
(tr-imag (imagpart tr))
109+
(known-vals nil)
110+
(unknown-vals nil))
111+
(loop :for (e1 e2) :on evals :by #'cddr
112+
:do (cond
113+
((and (realp e1) (realp e2))
114+
(assert (cl:= e1 e2))
115+
(push e1 known-vals)
116+
(decf tr-real e1))
117+
((and (complexp e1) (complexp e2))
118+
(assert (cl:= e1 (conjugate e2)))
119+
(push (complex (realpart e1)
120+
(abs (imagpart e1)))
121+
unknown-vals)
122+
(decf tr-real (realpart e1)))
123+
(t
124+
(error "unexpected eigenvalue pair"))))
125+
(format t "Re(tr) = ~A~%~
126+
Im(tr) = ~A~%~
127+
Known = ~A~%~
128+
UKnown = ~A~2%"
129+
tr-real
130+
tr-imag
131+
known-vals
132+
unknown-vals)
133+
(assert (< (abs tr-real) *junk-tol*))
134+
(format t "solving...~%")
135+
(loop :for sign :in (solve-plus-minus-sum
136+
(mapcar #'imagpart unknown-vals)
137+
tr-imag)
138+
:for unknown-val := (pop unknown-vals)
139+
:do (push (complex (realpart unknown-val)
140+
(* sign (imagpart unknown-val)))
141+
known-vals))
142+
143+
(format t "Known = ~A~%~
144+
UKnown = ~A~%~
145+
Tr = ~A~%~
146+
sum(eig) = ~A~%"
147+
known-vals
148+
unknown-vals
149+
tr
150+
(reduce #'+ known-vals))
151+
(assert (< (abs (- tr (reduce #'+ known-vals))) *junk-tol*))
152+
known-vals)))
153+
154+
155+
(defun solve-plus-minus-sum (a b)
156+
"Given a list of values A = (a1 a2 ... aN) and a value B, return a list of signs S = (s1 ... sN)---each of which is {-1, +1}---such that
157+
158+
S.A = B.
159+
"
160+
(labels ((rec (a s sum)
161+
(cond
162+
((null a)
163+
(cond
164+
((< (abs (- b sum)) *junk-tol*)
165+
(return-from solve-plus-minus-sum
166+
(reverse s)))
167+
((< (abs (+ b sum)) *junk-tol*)
168+
(return-from solve-plus-minus-sum
169+
(reverse (mapcar #'- s))))
170+
(t
171+
;; keep on truckin. return from REC and continue
172+
;; searching.
173+
)))
174+
(t
175+
(let ((ai (pop a)))
176+
(rec a (cons 1 s) (+ sum ai))
177+
(rec a (cons -1 s) (- sum ai)))))))
178+
(rec a nil 0)))

0 commit comments

Comments
 (0)