Commit 3467a214 authored by liam's avatar liam

Use letm instead of with-data, remove static objects in tests, add GSL

objects for histogram.


git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3283 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent ddee5998
;; Using GSL storage.
;; Liam Healy, Sun Mar 26 2006 - 16:32
;; Time-stamp: <2008-02-03 13:17:10EST data.lisp>
;; Time-stamp: <2008-02-03 23:34:26EST data.lisp>
;; $Id: $
(in-package :gsl)
......@@ -103,19 +103,42 @@
;;(LONG-DOUBLE . "_long_double")
(COMPLEX . "_complex")))
(defmacro data-go (type size-spec-is-list)
(defmacro data-go (type matrixp)
"Define the letm function for data types."
(let ((size-spec (if size-spec-is-list 'listp 'numberp))
(init-spec (if size-spec-is-list 'array-dimensions 'length)))
(if matrixp ; Matrix (two indices)
`(defgo ,type (size-or-initial &optional size-or-zero zero)
(if (and (numberp size-or-initial) (numberp size-or-zero))
;; Array dimensions are given as literal numbers for matrix
(list
`(make-data ',',type ,zero ,size-or-initial ,size-or-zero)
'free)
;; Determine at runtime whether first arg is initial or dimension
(let ((argsymb (gensym "ARG")))
(list
`(apply #'make-data ',',type nil
(if (and (numberp ,argsymb) (numberp ,size-or-zero))
(list ,argsymb ,size-or-zero)
(array-dimensions ,argsymb)))
'free
(lambda (symb)
`(unless (and (numberp ,argsymb) (numberp ,size-or-zero))
(setf (data ,symb) ,argsymb)))
(lambda () (list argsymb size-or-initial))))))
;; Vector or other one-index object
`(defgo ,type (size-or-initial &optional zero)
(if (numberp size-or-initial)
;; Size is given as literal number
(list
`(make-data ',',type ,zero ,size-or-initial)
'free)
;; Determine at runtime whether first arg is initial or size
(let ((argsymb (gensym "ARG")))
(list
`(make-data ',',type ,zero
(if (,',size-spec ,argsymb) ,argsymb (,',init-spec ,argsymb)))
(if (numberp ,argsymb) ,argsymb (length ,argsymb)))
'free
(lambda (symb)
`(unless (,',size-spec ,argsymb) (setf (data ,symb) ,argsymb)))
(lambda () (list argsymb size-or-initial)))))))
(lambda (symb) `(unless (numberp ,argsymb) (setf (data ,symb) ,argsymb)))
(lambda () (list argsymb size-or-initial))))))))
(defmacro defdata
(c-string cl-symbol cl-base-type
......
;; Permutations
;; Liam Healy, Sun Mar 26 2006 - 11:51
;; Time-stamp: <2008-01-21 11:12:01EST permutation.lisp>
;; Time-stamp: <2008-02-03 23:33:15EST permutation.lisp>
;; $Id: $
(in-package :gsl)
......@@ -69,8 +69,8 @@
(defun-gsl permutation-size (p)
"gsl_permutation_size"
(((pointer p) gsl-permutation-c))
:c-return :size ; FDL
:documentation
:c-return :size
:documentation ; FDL
"The size of the permutation p.")
(defun-gsl permutation-data (p)
......@@ -241,81 +241,68 @@
;;;; Examples and unit test
;;;;****************************************************************************
(defparameter *perm-1* (make-data 'permutation t 4))
(defparameter *perm-2* (make-data 'permutation t 4))
(defun generate-all-permutations (n)
"Generate all the permutations of n objects."
(letm ((perm (permutation n t)))
(loop collect (data perm 'list)
while (permutation-next perm))))
(lisp-unit:define-test permutation
(lisp-unit:assert-eql ;gsl-aref
2
(progn
(set-identity *perm-1*)
(gsl-aref *perm-1* 2)))
(letm ((perm-1 (permutation 4 t)))
(set-identity perm-1)
(gsl-aref perm-1 2)))
(lisp-unit:assert-equalp ;data
#(0 1 2 3)
(progn
(set-identity *perm-1*)
(data *perm-1*)))
(letm ((perm-1 (permutation 4 t)))
(set-identity perm-1)
(data perm-1)))
(lisp-unit:assert-equalp ;permutation-reverse
#(3 2 1 0)
(progn
(set-identity *perm-1*)
(data (permutation-reverse *perm-1*))))
(letm ((perm-1 (permutation 4 t)))
(set-identity perm-1)
(data (permutation-reverse perm-1))))
(lisp-unit:assert-equalp ;permutation-next, permutation-inverse
#(0 3 1 2)
(progn
(set-identity *perm-1*)
(permutation-next *perm-1*)
(permutation-next *perm-1*)
(permutation-next *perm-1*)
(permutation-inverse *perm-2* *perm-1*)
(data *perm-2*)))
(letm ((perm-1 (permutation 4 t)) (perm-2 (permutation 4 t)))
(set-identity perm-1)
(permutation-next perm-1)
(permutation-next perm-1)
(permutation-next perm-1)
(permutation-inverse perm-2 perm-1)
(data perm-2)))
(lisp-unit:assert-equalp ;swap-elements
#(0 3 2 1)
(progn
(set-identity *perm-1*)
(swap-elements *perm-1* 1 3)
(data *perm-1*)))
(letm ((perm-1 (permutation 4 t)))
(set-identity perm-1)
(swap-elements perm-1 1 3)
(data perm-1)))
(lisp-unit:assert-equalp ;permute-vector
#(33 44 11 22)
(progn
(set-identity *perm-1*)
(swap-elements *perm-1* 1 3)
(swap-elements *perm-1* 0 2)
(setf (data *intvec-1*) #(11 22 33 44))
(permute-vector *perm-1* *intvec-1*)
(data *intvec-1*)))
(letm ((perm-1 (permutation 4 t))
(intvec (vector-fixnum #(11 22 33 44))))
(set-identity perm-1)
(swap-elements perm-1 1 3)
(swap-elements perm-1 0 2)
(permute-vector perm-1 intvec)
(data intvec)))
(lisp-unit:assert-eql ;inversions
3
(progn
(set-identity *perm-1*)
(swap-elements *perm-1* 1 3)
(inversions *perm-1*)))
(letm ((perm-1 (permutation 4 t)))
(set-identity perm-1)
(swap-elements perm-1 1 3)
(inversions perm-1)))
(lisp-unit:assert-eql ;linear-cycles
3
(progn
(set-identity *perm-1*)
(swap-elements *perm-1* 1 3)
(linear-cycles *perm-1*)))
(letm ((perm-1 (permutation 4 t)))
(set-identity perm-1)
(swap-elements perm-1 1 3)
(linear-cycles perm-1)))
(lisp-unit:assert-eql ;canonical-cycles
2
(progn
(set-identity *perm-1*)
(swap-elements *perm-1* 1 3)
(swap-elements *perm-1* 0 2)
(canonical-cycles *perm-1*))))
;;;;****************************************************************************
;;;; Examples
;;;;****************************************************************************
#|
(with-data (perm permutation 4 t)
(loop collect (data perm 'list)
while (permutation-next perm)))
((0 1 2 3) (0 1 3 2) (0 2 1 3) (0 2 3 1) (0 3 1 2) (0 3 2 1) (1 0 2 3)
(1 0 3 2) (1 2 0 3) (1 2 3 0) (1 3 0 2) (1 3 2 0) (2 0 1 3) (2 0 3 1)
(2 1 0 3) (2 1 3 0) (2 3 0 1) (2 3 1 0) (3 0 1 2) (3 0 2 1) (3 1 0 2)
(3 1 2 0) (3 2 0 1) (3 2 1 0))
|#
(letm ((perm-1 (permutation 4 t)))
(set-identity perm-1)
(swap-elements perm-1 1 3)
(swap-elements perm-1 0 2)
(canonical-cycles perm-1))))
;; Eigenvectors and eigenvalues
;; Liam Healy, Sun May 21 2006 - 19:52
;; Time-stamp: <2008-02-03 13:22:47EST eigensystems.lisp>
;; Time-stamp: <2008-02-03 14:14:10EST eigensystems.lisp>
;; $Id: $
(in-package :gsl)
......@@ -182,21 +182,16 @@
;;;; Example
;;;;****************************************************************************
(defun eigenvalue-eigenvectors-example ()
(with-data (mat matrix-double (3 3))
(setf (data mat)
(letm ((evecs (matrix-double 3 3))
(evals (vector-double 3))
(ws (eigen-symmv 3))
(mat (matrix-double
#2A((20.0d0 -10.0d0 0.0d0)
(-10.0d0 30.0d0 0.0d0)
(0.0d0 0.0d0 40.0d0)))
(with-data (evecs matrix-double (3 3))
(with-data (evals vector-double 3)
(letm ((ws (eigen-symmv 3)))
(0.0d0 0.0d0 40.0d0)))))
(eigenvalues-eigenvectors-symmetric mat evals evecs ws)
(values (data evals) (data evecs)))))))
(values (data evals) (data evecs))))
#|
#(13.819660112501051d0 36.180339887498945d0 40.0d0)
......
;; Definition of GSLL system
;; Liam Healy
;; Time-stamp: <2008-01-28 22:24:15EST gsll.asd>
;; Time-stamp: <2008-02-03 22:30:09EST gsll.asd>
;; $Id: $
(asdf:defsystem "gsll"
......@@ -163,7 +163,8 @@
((:file "interpolation")
(:file "types")
(:file "lookup")
(:file "evaluation")))
(:file "evaluation")
(:file "spline-example")))
(:file "numerical-differentiation" :depends-on (init general))
(:file "chebyshev" :depends-on (init general))
(:file "series-acceleration" :depends-on (init))
......
;********************************************************
; file: structure.lisp
; description: The histogram structure.
; date: Mon Jan 1 2007 - 11:32
; author: Liam M. Healy
; modified: Tue Jan 2 2007 - 23:03
;********************************************************
;;; $Id: $
;; The histogram structure
;; Liam Healy, Mon Jan 1 2007 - 11:32
;; Time-stamp: <2008-02-03 21:53:32EST histogram.lisp>
;; $Id: $
(in-package :gsl)
;;; Define, make, copy histograms in one or two dimensions.
;;; Should this have a with-histogram macro, or be integrated into
;;; the data superclass (and so use with-data)?
(defgo histogram (size &rest from-or-ranges)
(list `(make-histogram ,size
,(when from-or-ranges
`(when (typep (first ,from-or-ranges) 'histogram)
(first ,from-or-ranges))))
'free
(when from-or-ranges
(lambda (sym)
`(unless (typep (first ,from-or-ranges) 'histogram)
(set-ranges ,sym ,@from-or-ranges))))))
(defclass histogram ()
((pointer
......@@ -26,7 +30,7 @@
"gsl_histogram_alloc"
(((number-of-bins object) :size))
:export nil
:index alloc
:index (letm histogram)
:c-return (cr :pointer)
:return ((assign-pointer object cr)))
......@@ -35,7 +39,7 @@
(((first (number-of-bins object)) :size)
((second (number-of-bins object)) :size))
:export nil
:index alloc
:index (letm histogram)
:c-return (cr :pointer)
:return ((assign-pointer object cr)))
......@@ -56,13 +60,17 @@
"gsl_histogram_free"
(((pointer object) :pointer))
:c-return (cr :pointer)
:return ((assign-pointer object cr)))
:return ((assign-pointer object cr))
:export nil
:index (letm histogram))
(defun-gsl free-histo-2 (object)
"gsl_histogram2d_free"
(((pointer object) :pointer))
:c-return (cr :pointer)
:return ((assign-pointer object cr)))
:return ((assign-pointer object cr))
:export nil
:index (letm histogram))
(defmethod free ((object histogram))
(histo-1d2d object free-histo))
......@@ -74,11 +82,11 @@
(((pointer histogram) :pointer)
((gsl-array ranges) :pointer) ((dim0 ranges) :size))
:export nil
:index set-ranges
:documentation
:index (letm histogram)
:documentation ; FDL
"Set the ranges of the existing histogram using
the gsl-vector of ranges. The values of the histogram
bins are reset to zero. The @code{ranges} array should contain the
bins are reset to zero. The ranges array should contain the
desired bin limits. The ranges can be arbitrary, subject to the
restriction that they are monotonically increasing.")
......@@ -88,7 +96,7 @@
((gsl-array x-ranges) :pointer) ((dim0 x-ranges) :size)
((gsl-array y-ranges) :pointer) ((dim0 y-ranges) :size))
:export nil
:index set-ranges)
:index (letm histogram))
(export 'set-ranges)
(defun set-ranges (histogram &rest ranges)
......@@ -109,16 +117,16 @@
(((pointer histogram) :pointer) (minimum :double) (maximum :double))
:export nil
:index set-ranges-uniform
:documentation
"Set the ranges of the existing histogram @var{h} to cover
the range @var{xmin} to @var{xmax} uniformly. The values of the
:documentation ; FDL
"Set the ranges of the existing histogram h to cover
the range xmin to xmax uniformly. The values of the
histogram bins are reset to zero. The bin ranges are shown in the table
below,
bin[0] corresponds to xmin <= x < xmin + d
bin[1] corresponds to xmin + d <= x < xmin + 2 d
......
bin[n-1] corresponds to xmin + (n-1)d <= x < xmax
where @math{d} is the bin spacing, @math{d = (xmax-xmin)/n}.")
where d is the bin spacing, d = (xmax-xmin)/n.")
;;; GSL documentation does not state what the return value for the
;;; C function means; assumed to be error code.
......@@ -133,15 +141,16 @@
(export 'set-ranges-uniform)
(defun set-ranges-uniform (histogram &rest limits)
"Set the ranges of the existing histogram @var{h} to cover
the range @var{xmin} to @var{xmax} uniformly. The values of the
;; FDL
"Set the ranges of the existing histogram h to cover
the range xmin to xmax uniformly. The values of the
histogram bins are reset to zero. The bin ranges are shown in the table
below,
bin[0] corresponds to xmin <= x < xmin + d
bin[1] corresponds to xmin + d <= x < xmin + 2 d
......
bin[n-1] corresponds to xmin + (n-1)d <= x < xmax
where @math{d} is the bin spacing, @math{d = (xmax-xmin)/n}."
where d is the bin spacing, d = (xmax-xmin)/n."
(histo-1d2d histogram set-ranges-uniform
((first limits) (second limits))
((first limits) (second limits)
......@@ -152,7 +161,7 @@
(((pointer destination) :pointer) ((pointer source) :pointer))
:export nil
:index copy
:documentation
:documentation ; FDL
"Copy the histogram source into the pre-existing
histogram destination, making the latter into
an exact copy of the former.
......@@ -172,7 +181,7 @@
(((pointer source) :pointer))
:export nil
:index clone
:documentation
:documentation ; FDL
"Create a new histogram which is an
exact copy of the histogram source, and return the pointer.")
......@@ -181,12 +190,15 @@
(((pointer source) :pointer))
:export nil
:index clone
:documentation
:documentation ; FDL
"Create a new histogram which is an
exact copy of the histogram source, and return the pointer.")
(export 'clone)
(defun clone (source)
;; FDL
"Create a new histogram which is an
exact copy of the histogram source, and return the pointer."
(histo-1d2d source clone))
(defun make-histogram (size &optional from)
......@@ -198,5 +210,3 @@
(unless (pointer ret)
(alloc ret))
ret))
......@@ -121,51 +121,61 @@
((first values) (second values))))
;;; Examples and unit test
(defparameter *sample-histogram*
(let ((histo (make-histogram 10)))
(set-ranges-uniform histo 0.0d0 10.0d0)
histo))
(defun setup-sample-histogram ()
(reset *sample-histogram*)
(increment *sample-histogram* 2.7d0)
(increment *sample-histogram* 6.9d0 2.0d0))
(lisp-unit:define-test histogram
(lisp-unit:assert-error
'gsl-warning
(progn
(setup-sample-histogram)
(increment *sample-histogram* -2.0d0)))
(letm ((histo (histogram 10)))
(set-ranges-uniform histo 0.0d0 10.0d0)
(increment histo -2.0d0)))
(lisp-unit:assert-first-fp-equal
"0.000000000000d+01"
(progn
(setup-sample-histogram)
(gsl-aref *sample-histogram* 1)))
(letm ((histo (histogram 10)))
(set-ranges-uniform histo 0.0d0 10.0d0)
(increment histo 2.7d0)
(increment histo 6.9d0 2.0d0)
(gsl-aref histo 1)))
(lisp-unit:assert-first-fp-equal
"0.100000000000d+01"
(progn
(setup-sample-histogram)
(gsl-aref *sample-histogram* 2)))
(letm ((histo (histogram 10)))
(set-ranges-uniform histo 0.0d0 10.0d0)
(increment histo 2.7d0)
(increment histo 6.9d0 2.0d0)
(gsl-aref histo 2)))
(lisp-unit:assert-first-fp-equal
"0.200000000000d+01"
(progn
(setup-sample-histogram)
(gsl-aref *sample-histogram* 6)))
(letm ((histo (histogram 10)))
(set-ranges-uniform histo 0.0d0 10.0d0)
(increment histo 2.7d0)
(increment histo 6.9d0 2.0d0)
(gsl-aref histo 6)))
(lisp-unit:assert-error
'gsl-error
(gsl-aref *sample-histogram* 16))
(lisp-unit:assert-first-fp-equal
"0.000000000000d+01"
(gsl-min-range *sample-histogram*))
(letm ((histo (histogram 10)))
(set-ranges-uniform histo 0.0d0 10.0d0)
(increment histo 2.7d0)
(increment histo 6.9d0 2.0d0)
(gsl-aref histo 16)))
(lisp-unit:assert-first-fp-equal
"0.100000000000d+02"
(gsl-max-range *sample-histogram*))
'("0.000000000000d+01" "0.100000000000d+02")
(lisp-unit:fp-values
(letm ((histo (histogram 10)))
(set-ranges-uniform histo 0.0d0 10.0d0)
(increment histo 2.7d0)
(increment histo 6.9d0 2.0d0)
(values (gsl-min-range histo) (gsl-max-range histo)))))
(lisp-unit:assert-eql
10
(bins *sample-histogram*))
(letm ((histo (histogram 10)))
(set-ranges-uniform histo 0.0d0 10.0d0)
(increment histo 2.7d0)
(increment histo 6.9d0 2.0d0)
(bins histo)))
(lisp-unit:assert-eql
5
(histogram-find *sample-histogram* 5.5d0)))
(letm ((histo (histogram 10)))
(set-ranges-uniform histo 0.0d0 10.0d0)
(increment histo 2.7d0)
(increment histo 6.9d0 2.0d0)
(histogram-find histo 5.5d0))))
;; Utility definitions
;; Liam Healy, Sun Dec 3 2006 - 10:21
;; Time-stamp: <2008-02-03 13:31:47EST utility.lisp>
;; Time-stamp: <2008-02-03 23:00:32EST utility.lisp>
;; $Id: $
(in-package :gsl)
......@@ -32,10 +32,11 @@
(let* (setters ; save setting forms
freers ; save freeing forms
assigns ; save global assignements
(*not-for-users* t)
(bnd ; rewrite the bindings with allocators
(mapcar
(lambda (b)
(*not-for-users* t))
;; Rewrite the bindings with allocators and push them on assigns
(dolist (b bindings)
(if (atom b)
(push b assigns)
(let ((defsymb (first b)) (defas (second b)))
(if (and (listp b) (listp defas)
(member (first defas) *letm-expanded-object*))
......@@ -49,22 +50,24 @@
;; Save the global assignments
(when assign (push (funcall assign) assigns))
;; Return the allocation form to be used in let
`(,defsymb ,allocator))
(push `(,defsymb ,allocator) assigns))
;; Not a GSL object, just define it with let
b)))
bindings)))
(push b assigns)))))
(if freers
;; There were GSL objects
`(let* (,@assigns ,@bnd)
`(let* (,@(reverse assigns))
;; put the body in an unwind-protect
(unwind-protect
(progn
,@setters ; setting the objects
,@(reverse setters) ; setting the objects
,@body)
,@freers)) ; freeing the objects
;; No GSL objects in the bindings, just make an ordinary let
`(let* (,@bindings) ,@body))))
;;; It might be nice to have a letm* or something where one go is
;;; dependent on another's being initialized.
(defun not-for-users (form)
(unless *not-for-users*
(error "The form ~a should be placed in the binding of a letm."
......
;; Interpolation allocation, initialization, and freeing.
;; Liam Healy, Sun Nov 4 2007 - 17:24
;; Time-stamp: <2008-02-03 13:17:11EST interpolation.lisp>
;; Time-stamp: <2008-02-03 22:55:28EST interpolation.lisp>
;; $Id: $
(in-package :gsl)
......@@ -14,7 +14,7 @@
'free-spline
(when ya
(lambda (symb)
`(initialize-interpolation ,symb ,xa-or-size ,ya (dim0 ,ya))))))
`(initialize-interpolation ,symb ,xa-or-size ,ya)))))
;;; Interpolation
(defun-gsl allocate-interpolation (type size)
......@@ -23,19 +23,19 @@
:c-return :pointer
:export nil
:index (letm interpolation)
:documentation
:documentation ; FDL
"Allocate an interpolation object of type for size data-points,
and return the pointer to it.")
(defun-gsl initialize-interpolation (interpolation xa ya size)
(defun-gsl initialize-interpolation (interpolation xa ya)
"gsl_interp_init"
((interpolation :pointer) (xa :pointer) (ya :pointer) (size :size))
:documentation
"Initialize the interpolation object @var{interp} for the
data (@var{xa},@var{ya}) where @var{xa} and @var{ya} are arrays of size
@var{size}. The interpolation object (@code{gsl_interp}) does not save
the data arrays @var{xa} and @var{ya} and only stores the static state
computed from the data. The @var{xa} data array is always assumed to be
((interpolation :pointer)
((gsl-array xa) :pointer) ((gsl-array ya) :pointer) ((dim0 xa) :size))
:documentation ; FDL
"Initialize the interpolation object interp for the
data (xa,ya) where xa and ya are gsl-arrays. The interpolation object does not save
the data arrays xa and ya and only stores the static state
computed from the data. The xa data array is always assumed to be
strictly ordered; the behavior for other arrangements is not defined.")
(defun-gsl free-interpolation (interpolation)
......@@ -53,31 +53,31 @@
'free-spline
(when ya
(lambda (symb)
`(initialize-spline ,symb ,xa-or-size ,ya (dim0 ,ya))))))
`(initialize-spline ,symb ,xa-or-size ,ya)))))
;;; Spline
(defun-gsl allocate-spline (type size)
"gsl_spline_alloc"
((type :pointer) (size :size))
:c-return :pointer
:documentation
:documentation ; FDL
"Allocate an interpolation object of type for size data-points,
and return the pointer to it.")
(defun-gsl initialize-spline (interpolation xa ya size)
(defun-gsl initialize-spline (interpolation xa ya)
"gsl_spline_init"
((interpolation :pointer) (xa :pointer) (ya :pointer) (size :size))
:documentation
"Initialize the interpolation object @var{interp} for the
data (@var{xa},@var{ya}) where @var{xa} and @var{ya} are arrays of size
@var{size}. The interpolation object (@code{gsl_spline}) does not save
the data arrays @var{xa} and @var{ya} and only stores the static state
computed from the data. The @var{xa} data array is always assumed to be
((interpolation :pointer)
((gsl-array xa) :pointer) ((gsl-array ya) :pointer) ((dim0 xa) :size))
:documentation ; FDL
"Initialize the interpolation object interp for the
data (xa,ya) where xa and ya are gsl-arrays. The spline object saves
the data arrays xa and ya computed from the data.
The xa data array is always assumed to be
strictly ordered; the behavior for other arrangements is not defined.")
(defun-gsl free-spline (interpolation)
"gsl_spline_free"
((interpolation :pointer))
:c-return :void
:documentation
"Frees the interpolation object @var{interp}.")
:documentation ; FDL
"Frees the spline object.")
;********************************************************
; file: spline-example.lisp
; description: Example spline
; date: Sat Nov 10 2007 - 21:18
; author: Liam Healy
; modified: Sat Nov 10 2007 - 22:10
;********************************************************
;;; $Id: $
;; Example spline
;; Liam Healy, Sat Nov 10 2007 - 21:18
;; Time-stamp: <2008-02-03 22:59:29EST spline-example.lisp>
;; $Id: $
(in-package :gsl)
;;; Possible future improvements:
;;; Direct double-float vectors.
;;; Use generic alloc/free defined in data.lisp
;;; Possible future improvement: direct double-float vectors.
(defun spline-example-arrays ()
(let ((xarr (make-array '(10) :element-type 'double-float))
......@@ -23,16 +17,11 @@
(defun spline-example ()
"The first example in Sec. 26.7 of the GSL manual."
(let ((acc (allocate-acceleration))
(spline (allocate-spline *cubic-spline-interpolation* 10)))
(multiple-value-bind (xarr yarr)
(spline-example-arrays)
(with-data (cxarr vector-double 10)
(with-data (cyarr vector-double 10)
(setf (data cxarr) xarr (data cyarr) yarr)
(initialize-spline spline (gsl-array cxarr) (gsl-array cyarr) 10)
(prog1
(letm ((acc (acceleration))
(cxarr (vector-double xarr))
(cyarr (vector-double yarr))
(spline (spline *cubic-spline-interpolation* cxarr cyarr)))
(loop for xi from (aref xarr 0) below (aref xarr 9) by 0.01d0
collect (list xi (evaluate-spline spline xi acc)))
(free-spline spline)
(free-acceleration acc)))))))
collect (list xi (evaluate-spline spline xi acc))))))
;********************************************************
; file: blas1.lisp
; description: BLAS level 1, Vector operations
; date: Wed Apr 26 2006 - 15:23
; author: Liam Healy
; modified: Mon Oct 8 2007 - 11:32
;********************************************************
;;; $Id: $
;; BLAS level 1, Vector operations
;; Liam Healy, Wed Apr 26 2006 - 15:23
;; Time-stamp: <2008-02-03 22:18:30EST blas1.lisp>
;; $Id: $
(in-package :gsl)
......@@ -22,13 +18,16 @@
(:documentation "Dot, or inner, product between vectors."))