;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; LUSH Lisp Universal Shell
;;;   Copyright (C) 2002 Leon Bottou, Yann Le Cun, AT&T Corp, NECI.
;;; Includes parts of TL3:
;;;   Copyright (C) 1987-1999 Leon Bottou and Neuristique.
;;; Includes selected parts of SN3.2:
;;;   Copyright (C) 1991-2001 AT&T Corp.
;;;
;;; This program is free software; you can redistribute it and/or modify
;;; it under the terms of the GNU General Public License as published by
;;; the Free Software Foundation; either version 2 of the License, or
;;; (at your option) any later version.
;;;
;;; This program is distributed in the hope that it will be useful,
;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;;; GNU General Public License for more details.
;;;
;;; You should have received a copy of the GNU General Public License
;;; along with this program; if not, write to the Free Software
;;; Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111, USA
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; $Id: cholesky.lsh,v 1.2 2003/06/04 18:40:44 leonb Exp $
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;; by Fu Jie Huang, May 2003
;; cholesky decomposition of real spd (symmetric positive definite) matrix

;;     [ 4  -2  4  2]
;; A = [-2  10 -2 -7]
;;     [ 4  -2  8  4]
;;     [ 2  -7  4  7]
;; A can decomposed as: A = U'*U
;;     [ 2  -1  2  1]
;; U = [     3  0 -2]
;;     [        2  1]
;;     [           1]

(libload "lapack/lapack-d")

(setq A [d [4 -2 4 2][-2 10 -2 -7][4 -2 8 4][2 -7 4 7]]) ;;initialized column-wise

(let ((uplo "U")
      (n [i@ 4])
      (lda [i@ 4])
      (info [i@ 0]))
  (dpotrf  uplo  n  A lda info))

(pretty A)

;; because of the row-column thing, you'll see U in A's lower triangle


