From 60c1f52d7f8804a93b1654978391e9886a19df8f Mon Sep 17 00:00:00 2001 From: "Liam M. Healy" <lhealy@common-lisp.net> Date: Sun, 30 Jan 2011 11:39:49 -0500 Subject: [PATCH] Remove invert-matrix to Antik, optional argument for LU-decomposition The user-friendly ("higher") functions layered on GSLL are now in Antik's math-high system, so the antik directory has been removed, and linear-algebra.lisp is enhanced and in Antik. #'LU-invert now takes 'inverse as an optional argument, with the sensible default. --- antik/linear-algebra.lisp | 94 --------------------------------------- gsll.asd | 8 +--- linear-algebra/lu.lisp | 9 ++-- 3 files changed, 8 insertions(+), 103 deletions(-) delete mode 100644 antik/linear-algebra.lisp diff --git a/antik/linear-algebra.lisp b/antik/linear-algebra.lisp deleted file mode 100644 index 6a6540f3..00000000 --- a/antik/linear-algebra.lisp +++ /dev/null @@ -1,94 +0,0 @@ -;; Define generic operators for GSL functions -;; Liam Healy 2011-01-01 09:51:47EST linear-algebra.lisp -;; Time-stamp: <2011-01-01 11:48:49EST linear-algebra.lisp> - -;; 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) - -(defmethod *i ((a grid:foreign-array) (b grid:foreign-array)) - (gsl:matrix-product a b)) - -(defmethod *i ((a grid:foreign-array) (b number)) - (gsl:scale b a)) - -(defmethod *i ((a number) (b grid:foreign-array)) - (gsl:scale a b)) - -(defmethod /i ((a grid:foreign-array) (b number)) - (gsl:scale (cl:/ b) a)) - -(defmethod +i ((a grid:foreign-array) (b grid:foreign-array)) - (gsl:elt+ a b)) - -(defmethod +i ((a grid:foreign-array) (b number)) - (gsl:elt+ a b)) - -(defmethod +i ((a number) (b grid:foreign-array)) - (gsl:elt+ b a)) - -(defmethod -i ((a grid:foreign-array) (b grid:foreign-array)) - (gsl:elt- a b)) - -(defmethod -i ((a grid:foreign-array) (b number)) - (gsl:elt+ a (cl:- b))) - -(defmethod -i ((a number) (b grid:foreign-array)) - (gsl:elt+ (- b) a)) - -(defmethod expt ((a grid:matrix) (b integer)) - (cond ((= 1 b) a) - ((plusp b) ; use addition-chain exponentiation - (* a (expt a (1- b)))) - ((minusp b) - (expt (invert-matrix a) (- b))) - ;; zero exponent; need to figure out non-square matrix - (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))))))) -|# - diff --git a/gsll.asd b/gsll.asd index 76a9c73f..0b755a7a 100644 --- a/gsll.asd +++ b/gsll.asd @@ -1,6 +1,6 @@ ;; Definition of GSLL system ;; Liam Healy -;; Time-stamp: <2011-01-30 10:00:04EST gsll.asd> +;; Time-stamp: <2011-01-30 10:41:57EST gsll.asd> ;; ;; Copyright 2006, 2007, 2008, 2009, 2010, 2011 Liam M. Healy ;; Distributed under the terms of the GNU General Public License @@ -266,11 +266,7 @@ :components ((cffi-grovel:grovel-file "mksa") (cffi-grovel:grovel-file "cgsm") - (:file export))) - (:module antik - :depends-on (init linear-algebra) - :components - ((:file "linear-algebra"))))) + (:file export))))) (asdf:defsystem-connection GSLL-tests :name "GSLL-tests" diff --git a/linear-algebra/lu.lisp b/linear-algebra/lu.lisp index 02b0a5e8..4ece7b4c 100644 --- a/linear-algebra/lu.lisp +++ b/linear-algebra/lu.lisp @@ -1,6 +1,6 @@ ;; LU decomposition ;; Liam Healy, Thu Apr 27 2006 - 12:42 -;; Time-stamp: <2011-01-12 00:44:43EST lu.lisp> +;; Time-stamp: <2011-01-30 11:00:48EST lu.lisp> ;; ;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy ;; Distributed under the terms of the GNU General Public License @@ -95,7 +95,10 @@ 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 grid:matrix) p inverse) +(defmfun LU-invert + ((LU grid:matrix) p + &optional (inverse + (grid:make-foreign-array element-type :dimensions (grid:dimensions LU)))) ("gsl_linalg" :complex "_LU_invert") (((mpointer LU) :pointer) ((mpointer p) :pointer) ((mpointer inverse) :pointer)) @@ -112,7 +115,7 @@ the same result more efficiently and reliably (consult any introductory textbook on numerical linear algebra for details).") - (defmfun LU-determinant ((LU grid:matrix) signum) +(defmfun LU-determinant ((LU grid:matrix) signum) ("gsl_linalg" :complex "_LU_det") (((mpointer LU) :pointer) (signum :int)) :definition :generic -- GitLab