(Arne Babenhauserheide)
2015-08-04: add example for cholesky decomposition. add example for cholesky decomposition.
diff --git a/examples/cholesky.w b/examples/cholesky.w new file mode 100644 --- /dev/null +++ b/examples/cholesky.w @@ -0,0 +1,50 @@ +#!/usr/bin/env sh +# -*- wisp -*- +exec guile -L $(dirname $(dirname $(realpath "$0"))) --language=wisp -e '(@@ (examples cholesky) main)' -s "$0" "$@" +; !# + +;; Cholesky decomposition, following https://de.wikipedia.org/wiki/Cholesky-Zerlegung#Pseudocode + +define-module : examples cholesky + . #:exports : cholesky! + +use-modules : guildhall ext foof-loop + +define : matrrix-ref X u v + list-ref (list-ref X u) v + +define : matrrix-set! X u v val + list-set! (list-ref X u) v val + + +define : cholesky! a + . "Modifies the square matirx a to contain the its cholesky decomposition. + +a is represented as list of lists." + let : : n : length a + loop + : for i : up-from 1 : to n + loop + : for j : up-from 1 : to i + let + : sum : matrix-ref a i j + when (>= j 1) + loop + : for k : up-from 1 : to (- j 1) + set! sum + - sum + * (matrix-ref a i k) (matrix-ref a j k) + cond + : > i j ; lower triangle + matrix-set! a i j + / sum : matrix-ref a j j + . a + : > sum 0 ; diagonal element + matrix-set! a i i : sqrt sum + . a + else + throw 'matrix-numerically-not-symmetric-positive-definite + + +define : main args + display : cholesky! '((1 2)(3 4))