Skip to content
Snippets Groups Projects
Commit f4bf4c1e authored by Liam M. Healy's avatar Liam M. Healy
Browse files

Move invert-matrix to antik/linear-algebra.lisp

parent 26271f71
No related branches found
No related tags found
No related merge requests found
;; Define generic operators for GSL functions ;; Define generic operators for GSL functions
;; Liam Healy 2011-01-01 09:51:47EST linear-algebra.lisp ;; Liam Healy 2011-01-01 09:51:47EST linear-algebra.lisp
;; Time-stamp: <2011-01-01 10:58:18EST linear-algebra.lisp> ;; Time-stamp: <2011-01-01 11:48:49EST linear-algebra.lisp>
;; $Id: $
;; Copyright 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 <http://www.gnu.org/licenses/>.
(in-package :antik) (in-package :antik)
...@@ -40,6 +55,40 @@ ...@@ -40,6 +55,40 @@
((plusp b) ; use addition-chain exponentiation ((plusp b) ; use addition-chain exponentiation
(* a (expt a (1- b)))) (* a (expt a (1- b))))
((minusp b) ((minusp b)
(expt (gsl:invert-matrix a) (- b))) (expt (invert-matrix a) (- b)))
;; zero exponent; need to figure out non-square matrix ;; zero exponent; need to figure out non-square matrix
(t (grid:identity-matrix (grid:dim0 a))))) (t (grid:identity-matrix (grid:dim0 a)))))
;;; Invert a matrix using LU
(export 'invert-matrix)
(defun invert-matrix (mat)
"Invert the matrix."
(let* ((dim (grid:dim0 mat))
(per (gsl:make-permutation dim))
(inv (grid:make-foreign-array 'double-float :dimensions (list dim dim))))
(gsl:lu-decomposition mat per)
(gsl:lu-invert mat per inv)))
#|
;;; Examples and unit test
;;; These are direct tests of matrix inversion
(save-test
lu
(grid:copy-to
(invert-matrix
(grid:make-foreign-array 'double-float
:dimensions '(2 2)
:initial-contents '(1.0d0 2.0d0 3.0d0 4.0d0)))))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST
#2A((-1.9999999999999998d0 1.0d0)
(1.4999999999999998d0 -0.49999999999999994d0)))
(MULTIPLE-VALUE-LIST
(GRID:COPY-TO
(INVERT-MATRIX
(GRID:MAKE-FOREIGN-ARRAY 'DOUBLE-FLOAT :DIMENSIONS '(2 2)
:INITIAL-CONTENTS
'((1.0d0 2.0d0) (3.0d0 4.0d0)))))))
|#
;; LU decomposition ;; LU decomposition
;; Liam Healy, Thu Apr 27 2006 - 12:42 ;; Liam Healy, Thu Apr 27 2006 - 12:42
;; Time-stamp: <2010-07-07 14:24:59EDT lu.lisp> ;; Time-stamp: <2011-01-01 11:47:45EST lu.lisp>
;; ;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy ;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License ;; Distributed under the terms of the GNU General Public License
;; ;;
;; This program is free software: you can redistribute it and/or modify ;; This program is free software: you can redistribute it and/or modify
...@@ -148,28 +148,6 @@ ...@@ -148,28 +148,6 @@
"Compute the sign or phase factor of the determinant of a matrix A, "Compute the sign or phase factor of the determinant of a matrix A,
det(A)/|det(A)|, from its LU decomposition, LU.") 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 (grid:make-foreign-array '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
(grid:copy-to
(invert-matrix
(grid:make-foreign-array 'double-float
:dimensions '(2 2)
:initial-contents '(1.0d0 2.0d0 3.0d0 4.0d0)))))
(generate-all-array-tests lu :doubles (generate-all-array-tests lu :doubles
(let ((matrix (array-default '(4 4))) (let ((matrix (array-default '(4 4)))
(vec (array-default '4 ))) (vec (array-default '4 )))
......
;; Regression test LU for GSLL, automatically generated ;; Regression test LU for GSLL, automatically generated
;; with some manual changes to the results ;; with some manual changes to the results
;; ;;
;; Copyright 2009 Liam M. Healy ;; Copyright 2009, 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License ;; Distributed under the terms of the GNU General Public License
;; ;;
;; This program is free software: you can redistribute it and/or modify ;; This program is free software: you can redistribute it and/or modify
...@@ -20,7 +20,7 @@ ...@@ -20,7 +20,7 @@
(in-package :gsl) (in-package :gsl)
;;; Answers (except for invert-matrix)inserted from linalg/test.c ;;; Answers inserted from linalg/test.c
;;; GSL has #define GSL_DBL_EPSILON 2.2204460492503131e-16 ;;; GSL has #define GSL_DBL_EPSILON 2.2204460492503131e-16
;;; which is 2x what double-float-epsilon is. ;;; which is 2x what double-float-epsilon is.
(LISP-UNIT:DEFINE-TEST LU (LISP-UNIT:DEFINE-TEST LU
...@@ -129,14 +129,4 @@ ...@@ -129,14 +129,4 @@
(MATRIX-PRODUCT-TRIANGULAR (MATRIX-PRODUCT-TRIANGULAR
MATRIX MATRIX
(MATRIX-PRODUCT-TRIANGULAR MATRIX X 1 :UPPER :NOTRANS :NONUNIT) (MATRIX-PRODUCT-TRIANGULAR MATRIX X 1 :UPPER :NOTRANS :NONUNIT)
1 :LOWER :NOTRANS :UNIT)))))))) 1 :LOWER :NOTRANS :UNIT)))))))))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST
#2A((-1.9999999999999998d0 1.0d0)
(1.4999999999999998d0 -0.49999999999999994d0)))
(MULTIPLE-VALUE-LIST
(GRID:COPY-TO
(INVERT-MATRIX
(GRID:MAKE-FOREIGN-ARRAY 'DOUBLE-FLOAT :DIMENSIONS '(2 2)
:INITIAL-CONTENTS
'((1.0d0 2.0d0) (3.0d0 4.0d0))))))))
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment