; ----------------------------------------------------------------------------
(define (cholesky A)
  (let* ((n (the fixnum (matrix-rows A)))
         (L (zeroes n)))
    (assume ((n fixnum) (k fixnum) (i fixnum) (j fixnum)
             (sum3 float) (init boolean)
             (A (vector fixnum fixnum (vector-of float)))
             (L (vector fixnum fixnum (vector-of float))))
    (do ((k 0 (+ k 1))) ((> k (- n 1)) #f)
        ;; First, calculate diagonal elements L_kk.
        (matrix-set! L k k
              (sqrt (- (matrix-ref A k k)
                       (let loop ((j 0) (sum3 0.0) (init #f))
                          (unless init
                             (set! init #t)
                             (set! sum3 (let ((x1 (matrix-ref L k j))) 
                                          (* x1 x1))))
                             (if (> j (- k 1))
                                 sum3
                                 (loop (+ j 1) (+ sum3 
                                                 (let ((x1 (matrix-ref L k (+ j 1))))
                                                    (* x1 x1))) 
                                                 init))))))
 
        ;; Then, all elements below a diagonal element, L_ik, i=k+1..n.
        (do ((i (+ k 1) (+ i 1)))
            ((> i (- n 1)) #f)
            (matrix-set! L i k
                  (/ (- (matrix-ref A i k)
                        (let loop ((j 0) (sum3 0.0) (init #f))
                          (unless init
                            (set! init #t)
                            (set! sum3 (* (matrix-ref L i j) (matrix-ref L k j))))
                             (if (> j (- k 1))
                                 sum3
                                 (loop (+ j 1) (+ sum3 (* (matrix-ref L i (+ j 1)) (matrix-ref L k (+ j 1)))) init))))
                     (matrix-ref L k k)))))
 
    ;; Return the calculated matrix L.
    L)))

; ----------------------------------------------------------------------------