Commit 60c1f52d authored by Liam M. Healy's avatar Liam M. Healy

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.
parent 476d01ee
;; 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)))))))
|#
;; 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"
......
;; 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
......
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