|
1 | 1 | ;;;; Thomas algorithm implementation in Common Lisp
|
2 | 2 |
|
| 3 | +(defmacro divf (place divisor) |
| 4 | + "Divides the value at place by divisor" |
| 5 | + `(setf ,place (/ ,place ,divisor))) |
| 6 | + |
| 7 | +(defun helper (v1 v2 v3 row) |
| 8 | + (- (svref v1 row) (* (svref v2 row) (svref v3 (1- row))))) |
| 9 | + |
3 | 10 | (defun thomas (diagonal-a diagonal-b diagonal-c last-column)
|
4 | 11 | "Returns the solutions to a tri-diagonal matrix non-destructively"
|
| 12 | + ;; We have to copy the inputs to ensure non-destructiveness |
5 | 13 | (let ((a (copy-seq diagonal-a))
|
6 | 14 | (b (copy-seq diagonal-b))
|
7 | 15 | (c (copy-seq diagonal-c))
|
8 | 16 | (d (copy-seq last-column)))
|
9 |
| - (setf (svref c 0) (/ (svref c 0) (svref b 0))) |
10 |
| - (setf (svref d 0) (/ (svref d 0) (svref b 0))) |
| 17 | + (divf (svref c 0) (svref b 0)) |
| 18 | + (divf (svref d 0) (svref b 0)) |
11 | 19 | (loop
|
12 | 20 | for i from 1 upto (1- (length a)) do
|
13 |
| - (setf |
14 |
| - (svref c i) |
15 |
| - (/ (svref c i) (- (svref b i) (* (svref a i) (svref c (1- i)))))) |
16 |
| - (setf |
17 |
| - (svref d i) |
18 |
| - (/ |
19 |
| - (- (svref d i) (* (svref a i) (svref d (1- i)))) |
20 |
| - (- (svref b i) (* (svref a i) (svref c (1- i))))))) |
| 21 | + (divf (svref c i) (helper b a c i)) |
| 22 | + (setf (svref d i) (/ (helper d a d i) (helper b a c i)))) |
21 | 23 | (loop
|
22 | 24 | for i from (- (length a) 2) downto 0 do
|
23 | 25 | (decf (svref d i) (* (svref c i) (svref d (1+ i)))))
|
|
28 | 30 | (defparameter diagonal-c #(4 5 0))
|
29 | 31 | (defparameter last-column #(7 5 3))
|
30 | 32 |
|
31 |
| -(format t "~{~f ~}" (coerce (thomas diagonal-a diagonal-b diagonal-c last-column) 'list)) |
| 33 | +;; should print 0.8666667 1.5333333 -0.26666668 |
| 34 | +(format t "~{~f ~}~%" (coerce (thomas diagonal-a diagonal-b diagonal-c last-column) 'list)) |
0 commit comments