;; Cholesky Decomposition
;; Liam Healy, Wed May 3 2006 - 16:38
;; Time-stamp: <2011-02-19 13:08:24EST cholesky.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
;;
;; 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 3 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, see .
(in-package :gsl)
;;; FDL
;;; A symmetric, positive definite square matrix A has a Cholesky
;;; decomposition into a product of a lower triangular matrix L and
;;; its transpose L^T,
;;; A = L L^T
;;; This is sometimes referred to as taking the square-root of a matrix. The
;;; Cholesky decomposition can only be carried out when all the eigenvalues
;;; of the matrix are positive. This decomposition can be used to convert
;;; the linear system A x = b into a pair of triangular systems
;;; (L y = b, L^T x = y), which can be solved by forward and
;;; back-substitution.
;;; GSL version 1.10 introduced functions for complex matrices.
(defmfun cholesky-decomposition ((A grid:matrix))
("gsl_linalg" :complex "_cholesky_decomp")
(((mpointer A) :pointer))
:definition :generic
:element-types :doubles
:inputs (A)
:outputs (A)
:documentation ; FDL
"Factorize the positive-definite square matrix A into the Cholesky
decomposition A = L L^T (real) or A = L L^H (complex). On output the
diagonal and lower triangular part of the input matrix A contain the
matrix L. The upper triangular part of the input matrix contains
L^T, the diagonal terms being identical for both L and L^T. If the
matrix is not positive-definite then the decomposition will fail,
returning the error input-domain.")
(defmfun cholesky-solve
((A grid:matrix) (b vector) &optional x-spec
&aux
(x (grid:ensure-foreign-array x-spec (grid:dimensions b) t)))
(("gsl_linalg" :complex "_cholesky_svx")
("gsl_linalg" :complex "_cholesky_solve"))
((((mpointer A) :pointer) ((mpointer b) :pointer))
(((mpointer A) :pointer) ((mpointer b) :pointer)
((mpointer x) :pointer)))
:definition :generic
:element-types :doubles
:inputs (A b)
:outputs (x)
:return ((or x b))
:documentation ; FDL
"Solve the system A x = b using the Cholesky
decomposition of A into the matrix given by
#'cholesky-decomposition. If x-spec is NIL (default), the solution
will replace b. If x-spec is T, then an array will be created and the
solution returned in it. If x-spec is a grid:foreign-array, the solution will
be returned in it.")
(defmfun cholesky-invert (cholesky)
"gsl_linalg_cholesky_invert"
(((mpointer cholesky) :pointer))
:inputs (cholesky)
:outputs (cholesky)
:gsl-version (1 12)
:documentation
"Compute the inverse of the matrix cholesky which must have been
previously computed by #'cholesky-decomposition. The inverse of the
original matrix is stored in cholesky on output.")
;;; Examples and unit test, from linalg/test.c
(defun test-cholesky-solve-dim (matrix)
"Solve the linear equation using Cholesky with the supplied matrix and
a right-hand side vector which is the reciprocal of one more than
the index."
(cholesky-solve
(cholesky-decomposition (copy matrix))
(create-rhs-vector (dim0 matrix))))
(defun test-cholesky-decomp-dim (matrix)
"Decompose using Cholesky and then multiply."
(let ((decomp (cholesky-decomposition (copy matrix))))
(dotimes (row (dim0 matrix) decomp)
(loop for col from (1+ row) below (dim1 matrix) do
(setf (grid:gref decomp row col) 0.0d0)))
(matrix-product decomp decomp nil 1.0d0 0.0d0 :notrans :trans)))
(defun test-cholesky-invert-dim (matrix)
"Invert using Cholesky decomposition"
(cholesky-invert (cholesky-decomposition (copy matrix))))
(save-test cholesky
(test-cholesky-solve-dim *hilb2*)
(test-cholesky-solve-dim *hilb3*)
(test-cholesky-solve-dim *hilb4*)
(test-cholesky-solve-dim *hilb12*)
(test-cholesky-decomp-dim *hilb2*)
(test-cholesky-decomp-dim *hilb3*)
(test-cholesky-decomp-dim *hilb4*)
(test-cholesky-decomp-dim *hilb12*)
(test-cholesky-invert-dim *hilb2*)
(test-cholesky-invert-dim *hilb3*)
(test-cholesky-invert-dim *hilb4*)
(test-cholesky-invert-dim *hilb12*))