2016-11-10: work on the cholesky decomposition

work on the cholesky decomposition

diff --git a/examples/cholesky.w b/examples/cholesky.w
--- a/examples/cholesky.w
+++ b/examples/cholesky.w
@@ -8,14 +8,29 @@ exec guile -L $(dirname $(dirname $(real
 define-module : examples cholesky
               . #:export : cholesky!
 
-use-modules : guildhall ext foof-loop
+use-modules : srfi srfi-42
 
-define : matrrix-ref X u v
-  list-ref (list-ref X u) v
+define : matrix-ref X row col
+  list-ref (list-ref X row) col
 
-define : matrrix-set! X u v val
-  list-set! (list-ref X u) v val
+define : matrix-set! X row col val
+  let : : sublist : list-ref X row
+    list-set! sublist col val
+    list-set! X row sublist
 
+define : matrix-transpose X
+   . "Swap columns and rows of a matrix"
+   list-ec (: outer (length X)) ; outer
+     list-ec (: inner (length (list-ref X outer))) ; inner
+       matrix-ref X inner outer
+
+define : matrix-multiply X Y
+       . "Calculate the matrix product of X and Y"
+       list-ec (: row (length Y))
+         list-ec (: col (length X))
+           sum-ec (: inner (length (list-ref Y row)))
+                   * : matrix-ref X inner col
+                       matrix-ref Y row inner
 
 define : cholesky! a
   . "Modifies the square matrix a to contain its cholesky decomposition.
@@ -24,23 +39,47 @@ sets a to g with a = ggT,
 
 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
+    do-ec (: i n)
+     begin
+      do-ec (: j (+ 1 i))
         let : : sum : matrix-ref a i j
+          ; format #t "n: ~A i: ~A j: ~A\n" n i j
           when : >= j 1
-            loop : : for k : up-from 1 : to {j - 1}
+            do-ec (: k {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
+                  matrix-ref a j j
             : > sum 0 ; diagonal element
               matrix-set! a i i : sqrt sum
-              . a
             else
-              throw 'matrix-numerically-not-symmetric-positive-definite
+              throw 'matrix-numerically-not-symmetric-positive-definite a
+      do-ec (: j (+ 1 i) n)
+        matrix-set! a i j 0
+    . a
 
 
 define : main args
-         display : cholesky! '((1 2)(2 4))
+   let : : X '((1 -1 1)
+               (-1 3 -.5)
+               (1 -.5 4))
+         format #t "cholesky\n"
+         display : cholesky! '((1 -1 1)
+                               (-1 3 -.5)
+                               (1 -.5 4))
+         newline
+         format #t "X\n"
+         display X
+         newline
+         cholesky! X
+         format #t "Lt·L\n"
+         display
+           matrix-multiply (matrix-transpose X) X
+         newline
+         display : matrix-transpose '((1 2)
+                                      (3 4))
+         newline
+         display : matrix-ref '((1 2)
+                                (3 4)) 0 1