Commit aa4aa130 authored by Liam Healy's avatar Liam Healy

LU decomposition ported to ffa

LU decomposition functions now support real and complex types, and the
example has been tested and works:
(invert-matrix
	     (make-array* '(2 2) 'double-float :initial-contents '(1.0d0 2.0d0 3.0d0 4.0d0)))
#2A((-1.9999999999999998d0 1.0d0)
    (1.4999999999999998d0 -0.49999999999999994d0))
parent 4a3fb0c5
;; Definition of GSLL system
;; Liam Healy
;; Time-stamp: <2008-08-10 18:27:45EDT gsll.asd>
;; Time-stamp: <2008-08-10 22:12:43EDT gsll.asd>
;; $Id$
(asdf:defsystem "gsll"
......@@ -83,7 +83,7 @@
(:file "blas2")
(:file "blas3" :depends-on (blas2))
(:file "exponential")
;(:file "lu")
(:file "lu")
;(:file "qr")
;(:file "qrpt")
;(:file "svd")
......
;; Macro for defining GSL functions.
;; Liam Healy 2008-04-16 20:49:50EDT defmfun.lisp
;; Time-stamp: <2008-08-10 18:00:49EDT defmfun.lisp>
;; Time-stamp: <2008-08-10 22:17:18EDT defmfun.lisp>
;; $Id$
(in-package :gsl)
......@@ -261,6 +261,7 @@
(:float *float-types*)
(:complex *complex-types*)
(:float-complex *float-complex-types*)
(:doubles *double-types*)
(t element-types)))))
(defun actual-gsl-function-name (base-name category type)
......@@ -273,16 +274,18 @@
(if (listp base-name)
(apply #'concatenate 'string
(substitute
(if (subtypep type 'complex) "u" "") :suffix
(if (subtypep type 'complex) "_complex" "") :complex
(substitute
(cl-gsl type (not blas) blas) :type
(if (subtypep type 'complex) "u" "") :suffix
(substitute
(if (subtypep type 'complex)
(cl-gsl (component-float-type type) nil blas)
"") :component-float-type
(cl-gsl type (not blas) blas) :type
(substitute
(string-downcase (symbol-name category)) :category
base-name))))))))
(if (subtypep type 'complex)
(cl-gsl (component-float-type type) nil blas)
"") :component-float-type
(substitute
(string-downcase (symbol-name category)) :category
base-name)))))))))
(defun number-class-from-type (type)
"Find a class name that contains this type."
......
;; Mapping of element type names
;; Liam Healy 2008-04-13 11:22:46EDT element-types.lisp
;; Time-stamp: <2008-07-20 22:34:41EDT element-types.lisp>
;; Time-stamp: <2008-08-10 22:17:19EDT element-types.lisp>
;; $Id$
;;; The different element type forms:
......@@ -299,6 +299,10 @@
(remove-if-not (lambda (tp) (subtypep tp 'complex)) *array-element-types*)
"All the supported complex array element types.")
(defparameter *double-types*
(list 'double-float '(complex double-float))
"All the supported double float element types.")
;;;;****************************************************************************
;;;; Types for CFFI (will eventually be in CFFI)
;;;;****************************************************************************
......
;; LU decomposition
;; Liam Healy, Thu Apr 27 2006 - 12:42
;; Time-stamp: <2008-03-30 18:51:01EDT lu.lisp>
;; Time-stamp: <2008-08-10 22:20:37EDT lu.lisp>
;; $Id$
(in-package :gsl)
;;; Not ported: those functions that use complex vectors or matrices
(defmfun LU-decomp (A p)
"gsl_linalg_LU_decomp"
(((pointer A) gsl-matrix-c) ((pointer p) gsl-permutation-c) (signum :int))
(defmfun LU-decomposition ((A matrix) p)
("gsl_linalg" :complex "_LU_decomp")
(((mpointer A) :pointer) ((mpointer p) :pointer) (signum :int))
:definition :generic
:inputs (A)
:outputs (A p)
:return (A)
:element-types :doubles
:documentation ; FDL
"Factorize the square matrix A into the LU
decomposition PA = LU. 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.
"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
p. The j-th column of the matrix P is given by the
......@@ -31,43 +32,59 @@
partial pivoting (Golub & Van Loan, Matrix Computations,
Algorithm 3.4.1).")
(defmfun LU-solve (LU p b x)
"gsl_linalg_LU_solve"
(((pointer LU) gsl-matrix-c) ((pointer p) gsl-permutation-c)
((pointer b) gsl-vector-c) ((pointer x) gsl-vector-c))
:invalidate (x)
(defmfun solve-LU ((LU matrix) p (b vector) (x vector))
("gsl_linalg" :complex "_LU_solve")
(((mpointer LU) :pointer) ((mpointer p) :pointer)
((mpointer b) :pointer) ((mpointer x) :pointer))
:definition :generic
:inputs (LU p b)
:outputs (x)
:return (x)
:element-types :doubles
:documentation ; FDL
"Solve the square system A x = b using the LU
decomposition of A into (LU, p) given by LU-decomp.")
(defmfun LU-svx (LU p x)
"gsl_linalg_LU_svx"
(((pointer LU) gsl-matrix-c) ((pointer p) gsl-permutation-c)
((pointer x) gsl-vector-c))
:invalidate (x)
(defmfun svx-LU ((LU matrix) p (x vector))
("gsl_linalg" :complex "_LU_svx")
(((mpointer LU) :pointer) ((mpointer p) :pointer)
((mpointer x) :pointer))
:definition :generic
:inputs (LU p)
:outputs (x)
:return (x)
:element-types :doubles
:documentation ; FLD
"Solve the square system A x = b in-place
using the LU decomposition of A into
(LU, p). On input x should contain the right-hand
side b, which is replaced by the solution on output.")
(defmfun LU-refine (A LU p b x residual)
"gsl_linalg_LU_refine"
(((pointer A) gsl-matrix-c) ((pointer LU) gsl-matrix-c)
((pointer p) gsl-permutation-c)
((pointer b) gsl-vector-c) ((pointer x) gsl-vector-c)
((pointer residual) gsl-vector-c))
:invalidate (x residual)
(defmfun refine-LU ((A matrix) LU p (b vector) (x vector) residual)
("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 p inverse)
"gsl_linalg_LU_invert"
(((pointer LU) gsl-matrix-c) ((pointer p) gsl-permutation-c)
((pointer inverse) gsl-matrix-c))
:invalidate (inverse)
(defmfun invert-LU ((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)
:return (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
......@@ -77,39 +94,39 @@
the same result more efficiently and reliably (consult any
introductory textbook on numerical linear algebra for details).")
(defgeneric LU-det (LU signum)
(:documentation ; FDL
"Compute the determinant of a matrix A from its LU
(defmfun determinant-LU ((LU matrix) signum)
("gsl_linalg" :complex "_LU_det")
(((mpointer LU) gsl-matrix-c) (signum :int))
:c-return :double
:definition :generic
:inputs (LU)
:element-types :doubles
: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-det ((LU matrix-double-float) signum)
"gsl_linalg_LU_det"
(((pointer LU) gsl-matrix-c) (signum :int))
:type :method
:c-return :double)
(defmfun LU-det ((LU matrix-complex) signum)
"gsl_linalg_complex_LU_det"
(((pointer LU) gsl-matrix-c) (signum :int))
:type :method
:c-return (ret gsl-complex)
:return ((complex-to-cl ret)))
(defmfun LU-lndet (LU)
"gsl_linalg_LU_lndet"
(((pointer LU) gsl-matrix-c))
diagonal elements of U and the sign of the row permutation signum.")
(defmfun log-determinant-LU ((LU matrix))
("gsl_linalg" :complex "_LU_lndet")
(((mpointer LU) gsl-matrix-c))
:c-return :double
:definition :generic
:inputs (LU)
:element-types :doubles
: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 signum)
"gsl_linalg_LU_sgndet"
(((pointer LU) gsl-matrix-c) (signum :int))
(defmfun LU-sgndet ((LU matrix) signum)
("gsl_linalg" :complex "_LU_sgndet")
(((mpointer LU) gsl-matrix-c) (signum :int))
:c-return :int
:definition :generic
:inputs (LU)
:element-types :doubles
:documentation ; FDL
"Compute the sign or phase factor of the determinant of a matrix A,
det(A)/|det(A)|, from its LU decomposition, LU.")
......@@ -121,20 +138,14 @@
(letm ((mmat (matrix-double-float mat))
(dim (array-dimension mat 0))
(per (permutation dim))
(inv (matrix-double-float dim dim)))
(lu-decomp mmat per)
(lu-invert mmat per inv)
(data inv)))
(inv (matrix-double-float (list dim dim))))
(LU-decomposition mmat per)
(invert-lu mmat per inv)
(cl-array inv)))
#|
(make-tests lu
(invert-matrix #2A((1.0d0 2.0d0) (3.0d0 4.0d0))))
(invert-matrix
(make-array* '(2 2) 'double-float :initial-contents '(1.0d0 2.0d0 3.0d0 4.0d0))))
|#
(LISP-UNIT:DEFINE-TEST LU
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST
#2A((-1.9999999999999998d0 1.0d0)
(1.4999999999999998d0 -0.49999999999999994d0)))
(MULTIPLE-VALUE-LIST
(INVERT-MATRIX #2A((1.0d0 2.0d0) (3.0d0 4.0d0))))))
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment