Thomas algorithm implementation in Common Lisp (#604)

Author: Trashtalk217

contents/thomas_algorithm/code/clisp/thomas.lisp

@@ -0,0 +1,34 @@
+;;;; Thomas algorithm implementation in Common Lisp
+
+(defmacro divf (place divisor)
+  "Divides the value at place by divisor"
+  `(setf ,place (/ ,place ,divisor)))
+
+(defun helper (v1 v2 v3 row)
+  (- (svref v1 row) (* (svref v2 row) (svref v3 (1- row)))))
+
+(defun thomas (diagonal-a diagonal-b diagonal-c last-column)
+  "Returns the solutions to a tri-diagonal matrix non-destructively"
+  ;; We have to copy the inputs to ensure non-destructiveness
+  (let ((a (copy-seq diagonal-a))
+         (b (copy-seq diagonal-b))
+         (c (copy-seq diagonal-c))
+         (d (copy-seq last-column)))
+    (divf (svref c 0) (svref b 0))
+    (divf (svref d 0) (svref b 0))
+    (loop
+      for i from 1 upto (1- (length a)) do
+      (divf (svref c i) (helper b a c i))
+      (setf (svref d i) (/ (helper d a d i) (helper b a c i))))
+    (loop
+      for i from (- (length a) 2) downto 0 do
+      (decf (svref d i) (* (svref c i) (svref d (1+ i)))))
+    d))
+
+(defparameter diagonal-a #(0 2 3))
+(defparameter diagonal-b #(1 3 6))
+(defparameter diagonal-c #(4 5 0))
+(defparameter last-column #(7 5 3))
+
+;; should print 0.8666667 1.5333333 -0.26666668
+(format t "~{~f ~}~%" (coerce (thomas diagonal-a diagonal-b diagonal-c last-column) 'list))

contents/thomas_algorithm/thomas_algorithm.md

@@ -129,6 +129,8 @@ You will find this algorithm implemented [in this project](https://scratch.mit.e
 [import, lang:"crystal"](code/crystal/thomas.cr)
 {% sample lang="kotlin" %}
 [import, lang:"kotlin"](code/kotlin/thomas.kt)
+{% sample lang="lisp" %}
+[import, lang:"lisp"](code/clisp/thomas.lisp)
 {% sample lang="ruby" %}
 [import, lang="ruby"](code/ruby/thomas.rb)
 {% sample lang="js" %}