#!/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)(2 4))