diff --git a/antik/linear-algebra.lisp b/antik/linear-algebra.lisp deleted file mode 100644 index 6a6540f3bdc76c8bfd26c36ad390ed5443af1ad5..0000000000000000000000000000000000000000 --- 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 76a9c73f0372cf49e6a327ec49351be6396270ed..0b755a7a2da7e9b3344cf358cb0c58438379d125 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 02b0a5e8b3a8352cc277a42ea34895a74a4469c4..4ece7b4c888a91afad005d7fdd8bb1586ffe3ed9 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