;; LU decomposition
;; Liam Healy, Thu Apr 27 2006 - 12:42
;; Time-stamp: <2010-02-05 09:15:11EST lu.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 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)
;;; /usr/include/gsl/gsl_linalg.h
(defmfun LU-decomposition
((A matrix) &optional (permutation (make-permutation (dim0 A))))
("gsl_linalg" :complex "_LU_decomp")
(((mpointer A) :pointer) ((mpointer permutation) :pointer)
(signum (:pointer :int)))
:definition :generic
:inputs (A)
:outputs (A permutation)
:return (A permutation signum)
:element-types :doubles
:documentation ; FDL
"Factorize the square matrix A into the LU decomposition PA = LU,
and return the sign of the permutation. On output the diagonal and
upper triangular part of the input matrix A contain the matrix U.
The lower triangular part of the input matrix (excluding the
diagonal) contains L. The diagonal elements of L are unity, and are
not stored.
The permutation matrix P is encoded in the permutation supplied as
the second argument and returned as the second value. The j-th
column of the matrix P is given by the k-th column of the identity
matrix, where k = p_j the j-th element of the permutation
vector. The sign of the permutation is returned as the second value;
it is the value (-1)^n, where n is the number of interchanges in the
permutation.
The algorithm used in the decomposition is Gaussian Elimination with
partial pivoting (Golub & Van Loan, Matrix Computations,
Algorithm 3.4.1).")
(defmfun LU-solve
((A matrix) (b vector) permutation &optional x-spec
&aux
(x (if (eq x-spec t)
(make-marray element-type :dimensions (dimensions b))
x-spec)))
(("gsl_linalg" :complex "_LU_svx")
("gsl_linalg" :complex "_LU_solve"))
((((mpointer A) :pointer)
((mpointer permutation) :pointer) ((mpointer b) :pointer))
(((mpointer A) :pointer)
((mpointer permutation) :pointer) ((mpointer b) :pointer)
((mpointer x) :pointer)))
:definition :generic
:element-types :doubles
:inputs (A b permutation)
:outputs (x b) ; depends on switch; both to be sure
:return ((or x b))
:documentation ; FDL
"Solve the square system A x = b using the LU
decomposition of A into (LU, p) given by LU-decomp.
If x-spec is nil, the solution will be computed in-place replacing b,
if it is T, an appropriate vector will be created and the solution
will be computed there. Otherwise it should be a supplied vector.")
(defmfun LU-refine
((A matrix) LU p (b vector) (x vector)
&optional (residual (make-marray element-type :dimensions (dim0 A))))
("gsl_linalg" :complex "_LU_refine")
(((mpointer A) :pointer) ((mpointer LU) :pointer)
((mpointer p) :pointer)
((mpointer b) :pointer) ((mpointer x) :pointer)
((mpointer residual) :pointer))
:definition :generic
:inputs (LU p)
:outputs (x residual)
:return (x)
:element-types :doubles
:documentation ; FDL
"Apply an iterative improvement to x, the solution of
A x = b, using the LU decomposition of A into (LU,p). The initial
residual r = A x - b is also computed and stored in residual. ")
(defmfun LU-invert ((LU matrix) p inverse)
("gsl_linalg" :complex "_LU_invert")
(((mpointer LU) :pointer) ((mpointer p) :pointer)
((mpointer inverse) :pointer))
:definition :generic
:inputs (LU p)
:outputs (inverse)
:element-types :doubles
:documentation ; FDL
"Compute the inverse of a matrix A from its LU
decomposition (LU,p), storing the result in the matrix inverse. The
inverse is computed by solving the system A x = b for each column of
the identity matrix. It is preferable to avoid direct use of the
inverse whenever possible, as the linear solver functions can obtain
the same result more efficiently and reliably (consult any
introductory textbook on numerical linear algebra for details).")
(defmfun LU-determinant ((LU matrix) signum)
("gsl_linalg" :complex "_LU_det")
(((mpointer LU) :pointer) (signum :int))
:definition :generic
:inputs (LU)
:element-types :doubles
:c-return :double
:documentation ; FDL
"Compute the determinant of a matrix from its LU
decomposition, LU. The determinant is computed as the product of the
diagonal elements of U and the sign of the row permutation signum.")
(defmfun LU-log-determinant ((LU matrix))
("gsl_linalg" :complex "_LU_lndet")
(((mpointer LU) :pointer))
:definition :generic
:inputs (LU)
:element-types :doubles
:c-return :double
:documentation ; FDL
"The logarithm of the absolute value of the
determinant of a matrix A, ln|det(A)|, from its LU decomposition,
LU. This function may be useful if the direct computation of the
determinant would overflow or underflow.")
(defmfun LU-sgndet ((LU matrix) signum)
("gsl_linalg" :complex "_LU_sgndet")
(((mpointer LU) :pointer) (signum (:pointer :int)))
:definition :generic
:inputs (LU)
:element-types :doubles
:c-return :int
:documentation ; FDL
"Compute the sign or phase factor of the determinant of a matrix A,
det(A)/|det(A)|, from its LU decomposition, LU.")
;;; Invert a matrix using LU
(export 'invert-matrix)
(defun invert-matrix (mat)
"Invert the matrix."
(let* ((dim (first (dimensions mat)))
(per (make-permutation dim))
(inv (make-marray 'double-float :dimensions (list dim dim))))
(LU-decomposition mat per)
(lu-invert mat per inv)))
;;; Examples and unit test
;;; These are direct tests of matrix inversion
(save-test
lu
(cl-array
(invert-matrix
(make-marray 'double-float
:dimensions '(2 2)
:initial-contents '(1.0d0 2.0d0 3.0d0 4.0d0)))))
(generate-all-array-tests lu :doubles
(let ((matrix (array-default '(4 4)))
(vec (array-default '4 )))
(multiple-value-bind (matrix perm)
(lu-decomposition matrix)
(let ((x (lu-solve matrix vec perm)))
(cl-array
(permute-inverse
perm
(matrix-product-triangular
matrix
(matrix-product-triangular matrix x 1 :upper :notrans :nonunit)
1 :lower :notrans :unit)))))))
;;; From linalg/test.c
(defun test-lu-solve-dim (matrix &optional vector)
"Solve the linear equation using LU with the supplied matrix and
a right-hand side vector which is the reciprocal of one more than
the index."
(let* ((dim (dim0 matrix))
(rhs (or vector
(create-rhs-vector dim (element-type matrix)))))
(multiple-value-bind (upper permutation signum)
(LU-decomposition (copy matrix))
(declare (ignore signum))
(let ((initial-solution (LU-solve upper rhs permutation T)))
(LU-refine matrix upper permutation rhs initial-solution)))))
(save-test lu
(test-lu-solve-dim *hilb2*)
(test-lu-solve-dim *hilb3*)
(test-lu-solve-dim *hilb4*)
(test-lu-solve-dim *hilb12*)
(test-lu-solve-dim *vander2*)
(test-lu-solve-dim *vander3*)
(test-lu-solve-dim *vander4*)
(test-lu-solve-dim *vander12*)
(test-lu-solve-dim (create-complex-matrix 7)))