Commit 2035b727 authored by liam's avatar liam

Introduction of #'with-gsl-objects macro and associated definitions,

and application to many chapters.  Most former with- macros have been
removed in favor of this macro.  All previous results have been
reproduced.  There is an experimental application to data (vectors,
matrices, etc.) but it is clumsy and setting doesn't work; this
requires more thought.


git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3276 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent 60f252ea
TODO
remove cffi-array
rename defun-gsl, defvariable.
Document with-gsl-objects
Features: Features:
- Designed so that usage of GSL is as if GSL were written in CL: - Designed so that usage of GSL is as if GSL were written in CL:
......
;; Chebyshev Approximations ;; Chebyshev Approximations
;; Liam Healy Sat Nov 17 2007 - 20:36 ;; Liam Healy Sat Nov 17 2007 - 20:36
;; Time-stamp: <2008-01-20 22:37:14EST chebyshev.lisp> ;; Time-stamp: <2008-01-26 18:17:32EST chebyshev.lisp>
;; $Id: $ ;; $Id: $
(in-package :gsl) (in-package :gsl)
;;; Needs with- macro to allocate, initialize, and free object.
;;;;**************************************************************************** ;;;;****************************************************************************
;;;; Creation and calculation of Chebyshev series ;;;; Creation and calculation of Chebyshev series
;;;;**************************************************************************** ;;;;****************************************************************************
(set-asf chebyshev allocate-chebyshev free-chebyshev initialize-chebyshev)
(defun-gsl allocate-chebyshev (order) (defun-gsl allocate-chebyshev (order)
"gsl_cheb_alloc" "gsl_cheb_alloc"
((order :size)) ((order :size))
:c-return :pointer :c-return :pointer
:export nil
:index (with-gsl-objects chebyshev)
:documentation :documentation
"Allocate a Chebyshev series of specified order "Allocate a Chebyshev series of specified order
and return a pointer to it.") and return a pointer to it.")
...@@ -23,6 +25,8 @@ ...@@ -23,6 +25,8 @@
"gsl_cheb_free" "gsl_cheb_free"
((chebyshev :pointer)) ((chebyshev :pointer))
:c-return :void :c-return :void
:export nil
:index (with-gsl-objects chebyshev)
:documentation :documentation
"Free a previously allocated Chebyshev series.") "Free a previously allocated Chebyshev series.")
...@@ -30,6 +34,8 @@ ...@@ -30,6 +34,8 @@
"gsl_cheb_init" "gsl_cheb_init"
((chebyshev :pointer) (function :pointer) ((chebyshev :pointer) (function :pointer)
(lower-limit :double) (upper-limit :double)) (lower-limit :double) (upper-limit :double))
:export nil
:index (with-gsl-objects chebyshev)
:documentation :documentation
"Compute the Chebyshev approximation for the function over the range "Compute the Chebyshev approximation for the function over the range
(lower-limit, upper-limit) to the previously specified order. The (lower-limit, upper-limit) to the previously specified order. The
...@@ -118,34 +124,28 @@ ...@@ -118,34 +124,28 @@
(defun-single chebyshev-step (x) (if (< x 0.5d0) 0.25d0 0.75d0)) (defun-single chebyshev-step (x) (if (< x 0.5d0) 0.25d0 0.75d0))
(defun chebyshev-table-example () (defun chebyshev-table-example ()
(let ((cheb (allocate-chebyshev 40)) (let ((steps 100))
(steps 100)) (with-gsl-objects ((chebyshev cheb 40 chebyshev-step 0.0d0 1.0d0))
(initialize-chebyshev cheb chebyshev-step 0.0d0 1.0d0) (dotimes (i steps)
(dotimes (i steps) (let ((x (coerce (/ i steps) 'double-float)))
(let ((x (coerce (/ i steps) 'double-float))) (format t "~&~a ~a ~a ~a"
(format t "~&~a ~a ~a ~a" x
x (chebyshev-step x)
(chebyshev-step x) (evaluate-chebyshev cheb x 10)
(evaluate-chebyshev cheb x 10) (evaluate-chebyshev cheb x)))))))
(evaluate-chebyshev cheb x))))
(free-chebyshev cheb)))
(defun chebyshev-point-example (x) (defun chebyshev-point-example (x)
(check-type x double-float) (check-type x double-float)
(let ((cheb (allocate-chebyshev 40)) (with-gsl-objects
(deriv (allocate-chebyshev 40)) ((chebyshev cheb 40 chebyshev-step 0.0d0 1.0d0)
(integ (allocate-chebyshev 40))) (chebyshev deriv 40)
(initialize-chebyshev cheb chebyshev-step 0.0d0 1.0d0) (chebyshev integ 40))
(derivative-chebyshev deriv cheb) (derivative-chebyshev deriv cheb)
(integral-chebyshev integ cheb) (integral-chebyshev integ cheb)
(prog1 (list
(list (evaluate-chebyshev cheb x)
(evaluate-chebyshev cheb x) (evaluate-chebyshev deriv x)
(evaluate-chebyshev deriv x) (evaluate-chebyshev integ x))))
(evaluate-chebyshev integ x))
(free-chebyshev cheb)
(free-chebyshev deriv)
(free-chebyshev integ))))
(lisp-unit:define-test chebyshev (lisp-unit:define-test chebyshev
(lisp-unit:assert-equal (lisp-unit:assert-equal
......
;; Using GSL storage. ;; Using GSL storage.
;; Liam Healy, Sun Mar 26 2006 - 16:32 ;; Liam Healy, Sun Mar 26 2006 - 16:32
;; Time-stamp: <2008-01-14 21:09:15 liam data.lisp> ;; Time-stamp: <2008-01-27 18:50:53EST data.lisp>
;; $Id: $ ;; $Id: $
(in-package :gsl) (in-package :gsl)
...@@ -198,6 +198,8 @@ ...@@ -198,6 +198,8 @@
(export '(make-data with-data)) (export '(make-data with-data))
(set-asf generic make-data free (setf data) 3)
(defgeneric alloc (object) (defgeneric alloc (object)
(:documentation "Allocate GSL data; used internally.")) (:documentation "Allocate GSL data; used internally."))
...@@ -210,13 +212,35 @@ ...@@ -210,13 +212,35 @@
(:method :after ((object gsl-data)) (setf (pointer object) nil))) (:method :after ((object gsl-data)) (setf (pointer object) nil)))
(defun make-data (type zero &rest size) (defun make-data (type zero &rest size)
"Make the GSL data object, including the allocation of space.
The user is responsible for calling #'free to free the foreign
memory when done."
(let ((obj
(make-instance
(make-symbol-from-strings *gsl-prefix* type)
:storage-size size)))
(if zero (calloc obj) (alloc obj))
obj))
(defmacro with-data ((symbol type size &optional zero) &body body)
"Allocate GSL data, bind to pointer,
and then deallocated it when done. If zero is T, zero the
contents when allocating."
`(let ((,symbol
(make-data ',type ,zero ,@(if (listp size) size (list size)))))
(unwind-protect
(progn ,@body)
(free ,symbol))))
#|
(defun make-data (type zero size)
"Make the GSL data object, including the allocation of space. "Make the GSL data object, including the allocation of space.
The user is responsible for calling #'free to free the foreign The user is responsible for calling #'free to free the foreign
memory when done." memory when done."
(let ((obj (let ((obj
(make-instance (make-instance
(make-symbol-from-strings *gsl-prefix* type) (make-symbol-from-strings *gsl-prefix* type)
:storage-size size))) :storage-size (if (listp size) size (list size)))))
(if zero (calloc obj) (alloc obj)) (if zero (calloc obj) (alloc obj))
obj)) obj))
...@@ -225,10 +249,11 @@ ...@@ -225,10 +249,11 @@
and then deallocated it when done. If zero is T, zero the and then deallocated it when done. If zero is T, zero the
contents when allocating." contents when allocating."
`(let ((,symbol `(let ((,symbol
(make-data ',type ,zero ,@(if (listp size) size (list size))))) (make-data ',type ,zero ,size)))
(unwind-protect (unwind-protect
(progn ,@body) (progn ,@body)
(free ,symbol)))) (free ,symbol))))
|#
;;;;**************************************************************************** ;;;;****************************************************************************
;;;; Getting values into CL ;;;; Getting values into CL
......
;; Vectors ;; Vectors
;; Liam Healy, Sun Mar 26 2006 - 11:51 ;; Liam Healy, Sun Mar 26 2006 - 11:51
;; Time-stamp: <2008-01-15 19:01:37 liam vector.lisp> ;; Time-stamp: <2008-01-26 22:20:43EST vector.lisp>
;; $Id: $ ;; $Id: $
(in-package :gsl) (in-package :gsl)
...@@ -469,6 +469,14 @@ same length.") ...@@ -469,6 +469,14 @@ same length.")
(setf (data vec) #(-3.21d0 1.0d0 12.8d0)) (setf (data vec) #(-3.21d0 1.0d0 12.8d0))
(data vec)) (data vec))
(with-gsl-objects ((generic vec 'vector-double nil 3))
(setf (data vec) #(-3.21d0 1.0d0 12.8d0))
(data vec))
(with-data (vec vector 3) (set-basis vec 1) (with-data (vec vector 3) (set-basis vec 1)
(format t "~&~a ~a ~a" (format t "~&~a ~a ~a"
(gsl-aref vec 0) (gsl-aref vec 1) (gsl-aref vec 2))) (gsl-aref vec 0) (gsl-aref vec 1) (gsl-aref vec 2)))
......
;; Eigenvectors and eigenvalues ;; Eigenvectors and eigenvalues
;; Liam Healy, Sun May 21 2006 - 19:52 ;; Liam Healy, Sun May 21 2006 - 19:52
;; Time-stamp: <2008-01-21 11:30:57EST eigensystems.lisp> ;; Time-stamp: <2008-01-27 18:53:17EST eigensystems.lisp>
;; $Id: $ ;; $Id: $
(in-package :gsl) (in-package :gsl)
;;; Eventually make a with-* macro for the eigen-herm functions.
;;;;****************************************************************************
;;;; Workspace macro
;;;;****************************************************************************
(export 'with-eigensystem-workspace)
(defmacro with-eigensystem-workspace
((ws n &optional vectors hermitian) &body body)
"Create a workspace to use in finding eigenvalues and eigenvectors.
If eigenvectors are desired, set 'vectors to T."
`(let ((,ws ,(if vectors
(if hermitian
`(eigen-hermv-alloc ,n)
`(eigen-symmv-alloc ,n))
(if hermitian
`(eigen-herm-alloc ,n)
`(eigen-symm-alloc ,n)))))
(check-null-pointer ,ws :ENOMEM "workspace allocation failed")
(unwind-protect
(progn ,@body)
,(if vectors
(if hermitian
`(eigen-hermv-free ,n)
`(eigen-symmv-free ,ws))
(if hermitian
`(eigen-herm-free ,n)
`(eigen-symm-free ,ws))))))
;;;;**************************************************************************** ;;;;****************************************************************************
;;;; Real Symmetric Matrices ;;;; Real Symmetric Matrices
;;;;**************************************************************************** ;;;;****************************************************************************
(set-asf eigen-symm eigen-symm-alloc eigen-symm-free)
(defun-gsl eigen-symm-alloc (n) (defun-gsl eigen-symm-alloc (n)
"gsl_eigen_symm_alloc" ((n :size)) "gsl_eigen_symm_alloc" ((n :size))
:c-return :pointer
:export nil
:index '(with-gsl-objects eigen-symm)
:documentation ; FDL :documentation ; FDL
"Allocate a workspace for computing eigenvalues of "Allocate a workspace for computing eigenvalues of
n-by-n real symmetric matrices. The size of the workspace n-by-n real symmetric matrices. The size of the workspace
is O(2n)." is O(2n).")
:c-return :pointer)
(defun-gsl eigen-symm-free (w) (defun-gsl eigen-symm-free (w)
"gsl_eigen_symm_free" ((w :pointer)) "gsl_eigen_symm_free" ((w :pointer))
:c-return :void :c-return :void
:export nil
:index '(with-gsl-objects eigen-symm)
:documentation ; FDL :documentation ; FDL
"Free the memory associated with the workspace w.") "Free the memory associated with the workspace w.")
(set-asf eigen-symmv eigen-symmv-alloc eigen-symmv-free)
(defun-gsl eigen-symmv-alloc (n) (defun-gsl eigen-symmv-alloc (n)
"gsl_eigen_symmv_alloc" ((n :size)) "gsl_eigen_symmv_alloc" ((n :size))
:c-return :pointer
:export nil
:index '(with-gsl-objects eigen-symmv)
:documentation ; FDL :documentation ; FDL
"Allocate a workspace for computing eigenvalues and "Allocate a workspace for computing eigenvalues and
eigenvectors of n-by-n real symmetric matrices. The size of eigenvectors of n-by-n real symmetric matrices. The size of
the workspace is O(4n)." the workspace is O(4n).")
:c-return :pointer)
(defun-gsl eigen-symmv-free (w) (defun-gsl eigen-symmv-free (w)
"gsl_eigen_symmv_free" ((w :pointer)) "gsl_eigen_symmv_free" ((w :pointer))
:index '(with-gsl-objects eigen-symmv)
:c-return :void :c-return :void
:documentation ; FDL :documentation ; FDL
"Free the memory associated with the workspace w.") "Free the memory associated with the workspace w.")
...@@ -81,7 +62,8 @@ ...@@ -81,7 +62,8 @@
(defun-gsl eigenvalues-eigenvectors-symmetric (A eval evec ws) (defun-gsl eigenvalues-eigenvectors-symmetric (A eval evec ws)
"gsl_eigen_symmv" "gsl_eigen_symmv"
((A gsl-matrix-c) (eval gsl-vector-c) (evec gsl-matrix-c) (ws :pointer)) (((pointer A) :pointer) ((pointer eval) :pointer)
((pointer evec) :pointer) (ws :pointer))
:documentation ; FDL :documentation ; FDL
"The eigenvalues and eigenvectors of the real "The eigenvalues and eigenvectors of the real
symmetric matrix A. Additional workspace of the appropriate size symmetric matrix A. Additional workspace of the appropriate size
...@@ -100,6 +82,8 @@ ...@@ -100,6 +82,8 @@
;;;; Complex Hermitian Matrices ;;;; Complex Hermitian Matrices
;;;;**************************************************************************** ;;;;****************************************************************************
(set-asf eigen-herm eigen-herm-alloc eigen-herm-free)
(defun-gsl eigen-herm-alloc (n) (defun-gsl eigen-herm-alloc (n)
"gsl_eigen_herm_alloc" ((n :size)) "gsl_eigen_herm_alloc" ((n :size))
:documentation ; FDL :documentation ; FDL
...@@ -128,6 +112,8 @@ ...@@ -128,6 +112,8 @@
:invalidate (eval A) :invalidate (eval A)
:return (eval)) :return (eval))
(set-asf eigen-hermv eigen-hermv-alloc eigen-hermv-free)
(defun-gsl eigen-hermv-alloc (n) (defun-gsl eigen-hermv-alloc (n)
"gsl_eigen_hermv_alloc" ((n :size)) "gsl_eigen_hermv_alloc" ((n :size))
:documentation ; FDL :documentation ; FDL
...@@ -196,19 +182,23 @@ ...@@ -196,19 +182,23 @@
;;;; Example ;;;; Example
;;;;**************************************************************************** ;;;;****************************************************************************
#|
(with-data (mat matrix (3 3)) (defun eigenvalue-eigenvectors-example ()
(setf (data mat) (with-data (mat matrix-double (3 3))
#2A((20.0d0 -10.0d0 0.0d0) (setf (data mat)
(-10.0d0 30.0d0 0.0d0) #2A((20.0d0 -10.0d0 0.0d0)
(0.0d0 0.0d0 40.0d0))) (-10.0d0 30.0d0 0.0d0)
(with-data (evecs matrix (3 3)) (0.0d0 0.0d0 40.0d0)))
(with-data (evals vector 3) (with-data (evecs matrix-double (3 3))
(with-eigensystem-workspace (ws 3 t) (with-data (evals vector-double 3)
(eigenvalues-eigenvectors-symmetric mat evals evecs ws) (with-gsl-objects ((eigen-symmv ws 3))
(values (data evals) (data evecs)))))) (eigenvalues-eigenvectors-symmetric mat evals evecs ws)
(values (data evals) (data evecs)))))))
#|
#(13.819660112501051d0 36.180339887498945d0 40.0d0) #(13.819660112501051d0 36.180339887498945d0 40.0d0)
#2A((0.8506508083520399d0 -0.5257311121191336d0 0.0d0) #2A((0.8506508083520399d0 -0.5257311121191336d0 0.0d0)
(0.5257311121191336d0 0.8506508083520399d0 0.0d0) (0.5257311121191336d0 0.8506508083520399d0 0.0d0)
......
;******************************************************** ;; Discrete Hankel Transforms.
; file: hankel.lisp ;; Liam Healy, Sat Dec 8 2007 - 16:50
; description: Discrete Hankel Transforms. ;; Time-stamp: <2008-01-26 18:36:00EST hankel.lisp>
; date: Sat Dec 8 2007 - 16:50 ;; $Id: $
; author: Liam Healy
; modified: Sat Dec 8 2007 - 18:38
;********************************************************
;;; $Id: $
(in-package :gsl) (in-package :gsl)
;;; Everything compiles, but not tested -- need example. ;;; Everything compiles, but not tested -- need example.
(set-asf hankel allocate-hankel free-hankel init-hankel)
(defun-gsl allocate-hankel (size) (defun-gsl allocate-hankel (size)
"gsl_dht_alloc" "gsl_dht_alloc"
((size :size)) ((size :size))
:c-return :pointer :c-return :pointer
:export nil
:index (with-gsl-objects hankel)
:documentation :documentation
"Allocate a Discrete Hankel transform object of size @var{size}.") "Allocate a Discrete Hankel transform object of size @var{size}.")
(defun-gsl init-hankel (hankel nu xmax) (defun-gsl init-hankel (hankel nu xmax)
"gsl_dht_new" "gsl_dht_new"
((hankel :pointer) (nu :double) (xmax :double)) ((hankel :pointer) (nu :double) (xmax :double))
:export nil
:index (with-gsl-objects hankel)
:documentation :documentation
"Initialize the transform @var{t} for the given values of "Initialize the transform @var{t} for the given values of
@var{nu} and @var{x}.") @var{nu} and @var{x}.")
;; This seems redundant; does it effectively combine allocate-hankel
;; and init-hankel? In this case we don't really need, and shouldn't
;; export it.
(defun-gsl new-hankel (size nu xmax) (defun-gsl new-hankel (size nu xmax)
"gsl_dht_new" "gsl_dht_new"
((size :size) (nu :double) (xmax :double)) ((size :size) (nu :double) (xmax :double))
...@@ -38,6 +43,8 @@ ...@@ -38,6 +43,8 @@
"gsl_dht_free" "gsl_dht_free"
((hankel :pointer)) ((hankel :pointer))
:c-return :void :c-return :void
:export nil
:index (with-gsl-objects hankel)
:documentation :documentation
"Free the Hankel object.") "Free the Hankel object.")
......
...@@ -206,7 +206,7 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre>.</p> ...@@ -206,7 +206,7 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre>.</p>
</tr> </tr>
<tr> <tr>
<td><a href="Vectors-and-Matrices.html">Vectors and Matrices</a></td> <td><a href="Vectors-and-Matrices.html">Vectors and Matrices</a></td>
<td>Done for integer, real, double, complex</td> <td>Done for integer, real, double, complex; subvectors don't work</td>
</tr> </tr>
<tr> <tr>
<td><a href="Permutations.html">Permutations</a></td> <td><a href="Permutations.html">Permutations</a></td>
...@@ -374,7 +374,7 @@ and expect a few failures on amd64 (four), none? on i386. Often, the ...@@ -374,7 +374,7 @@ and expect a few failures on amd64 (four), none? on i386. Often, the
<!-- Created: Feb 25 2005 --> <!-- Created: Feb 25 2005 -->
<!-- hhmts start --> <!-- hhmts start -->
<small> <small>
Time-stamp: <2008-01-21 17:58:20EST index.html> Time-stamp: <2008-01-26 22:03:36EST index.html>
</small> </small>
<!-- hhmts end --> <!-- hhmts end -->
</div> </div>
......
;******************************************************** ;; Utility definitions
; file: utility.lisp ;; Liam Healy, Sun Dec 3 2006 - 10:21
; description: Utility definitions ;; Time-stamp: <2008-01-27 13:28:29EST utility.lisp>
; date: Sun Dec 3 2006 - 10:21 ;; $Id: $
; author: Liam M. Healy
; modified: Sun Dec 3 2006 - 10:22
;********************************************************
;;; $Id: $
(in-package :gsl) (in-package :gsl)
...@@ -14,3 +10,87 @@ ...@@ -14,3 +10,87 @@
(intern (format nil "~{~:@(~a~)~^-~}" (mapcar #'string strings)))) (intern (format nil "~{~:@(~a~)~^-~}" (mapcar #'string strings))))
(defparameter *gsl-prefix* 'gsl) (defparameter *gsl-prefix* 'gsl)
(defparameter *gsl-objects* nil)
(defstruct gsl-object name allocate free set num-alloc-args)
(defmacro set-asf (name allocate free &optional set (num-alloc-args 1))
"Define an expansion of this name in #'with-gsl-objects
and use the appropriate functions for allocating, freeing, and setting."
`(eval-when (:compile-toplevel :load-toplevel :execute)
(pushnew
(make-gsl-object :name ',name :allocate ',allocate :free ',free
:set ',set :num-alloc-args ,num-alloc-args)
*gsl-objects*
:test #'equalp)))
;; (set-asf poly-complex complex-workspace-alloc complex-workspace-free nil 1)
;; (set-asf mfminimizer allocate-mfminimizer free-mfminimizer set-mfminimizer 2)
(defun gslo-letform (specs specfns)
"Create the let form for GSL objects."
(mapcar
(lambda (s sf)
`(,(second s)
(,(gsl-object-allocate sf)
,@(subseq (cddr s) 0 (gsl-object-num-alloc-args sf)))))
specs specfns))
(defun gslo-setform (specs specfns)
"Create the set forms for GSL objects."
(loop for s in specs for sf in specfns
for addl-set-args
= (subseq (cddr s) (gsl-object-num-alloc-args sf))
when (and (gsl-object-set sf) addl-set-args)
collect
`(,(gsl-object-set sf)
;; presume the thing is the first arg to the set function,
;; then add the other arguments
,@(cons (second s) addl-set-args))))
(export 'with-gsl-objects)
(defmacro with-gsl-objects (specs &body body)
"Allocate, optionally set, and free any foreign object
used in GSL calculations."
(let ((specfns
(mapcar
(lambda (s)
(or (find (first s) *gsl-objects* :key #'gsl-object-name)
(error "Could not find GSL object named ~s" (first s))))
specs)))
`(let (,@(gslo-letform specs specfns))
(unwind-protect
(progn ; call set functions
,@(gslo-setform specs specfns)
,@body)
,@(loop for s in specs for sf in specfns
collect `(,(gsl-object-free sf) ,(second s)))))))
;;; roots-multi needs lambda list
;;; (solver solver-type function initial)
;;; allocate: solver-type (dim0 ,initial)
;;; set: solver function initial
#|
(vector-double vec init)
->
(make-data vector-double nil (length init))
(setf (data vec) init)
;;; Example
(set-asf chebyshev allocate-chebyshev free-chebyshev initialize-chebyshev)
(with-gsl-objects ((chebyshev cheb 40 chebyshev-step 0.0d0 1.0d0))
...)
(LET ((CHEB (ALLOCATE-CHEBYSHEV 40)))
(UNWIND-PROTECT
(PROGN
(INITIALIZE-CHEBYSHEV CHEB CHEBYSHEV-STEP 0.0d0 1.0d0)
...)
(FREE-CHEBYSHEV CHEB)))
|#