Commit 520e168e authored by Sumant Oemrawsingh's avatar Sumant Oemrawsingh
Browse files

Merge branch 'master' of ssh://repo.or.cz/srv/git/gsll

parents b966b76d 360080b1
;; Permutations
;; Liam Healy, Sun Mar 26 2006 - 11:51
;; Time-stamp: <2009-06-06 09:49:18EDT permutation.lisp>
;; Time-stamp: <2009-11-14 10:52:43EST permutation.lisp>
;; $Id$
(in-package :gsl)
......@@ -106,6 +106,7 @@
"A pointer to the array of elements in the
permutation p.")
(export 'validp)
(defgeneric validp (object)
(:documentation ; FDL
"Check that the object p is valid."))
......
;; Wavetable and workspace object definitions.
;; Sumant Oemrawsingh, Sun Oct 25 2009 - 16:35
;; Time-stamp: <2009-11-07 11:24:02EST wavetable-workspace.lisp>
;; Time-stamp: <2009-11-11 19:01:48EST wavetable-workspace.lisp>
(in-package :gsl)
......@@ -212,8 +212,8 @@
(make-fft-real-wavetable-single-float dimension)))
((eql element-type 'double-float)
(if half-complex
(make-fft-real-wavetable-double-float dimension)
(make-fft-half-complex-wavetable-double-float dimension)))
(make-fft-half-complex-wavetable-double-float dimension)
(make-fft-real-wavetable-double-float dimension)))
((equal element-type '(complex single-float))
(make-fft-complex-wavetable-single-float dimension))
((equal element-type '(complex double-float))
......
;; Expand the body of a defmfun
;; Liam Healy 2009-04-13 22:07:13EDT body-expand.lisp
;; Time-stamp: <2009-05-03 22:18:13EDT body-expand.lisp>
;; Time-stamp: <2009-11-14 10:19:46EST body-expand.lisp>
;; $Id: $
(in-package :gsl)
......@@ -27,19 +27,6 @@
(make-symbol (string (st-symbol sd)))
(st-type sd))))))
(defun substitute-symbols (substlist st-list)
"Substitute in the new symbols for the old in the st list."
(loop for pair in substlist
with answer = st-list
do
(setf answer
(subst
(st-symbol (second pair))
(first pair)
answer
:test 'st-symbol))
finally (return answer)))
;;; This function should never be called even when FSBV is absent,
;;; because the potential callers should all have #-fsbv
;;; conditionalization. It is here just so that body-expand can
......
;; Define examples.
;; Liam Healy 2008-09-07 21:00:48EDT generate-tests.lisp
;; Time-stamp: <2009-02-11 22:44:32EST generate-examples.lisp>
;; Time-stamp: <2009-11-14 10:22:05EST generate-examples.lisp>
;; $Id: $
;;; Define examples that can be displayed by users with the function
......@@ -75,15 +75,6 @@
-82 75 -30 -71 44 48 88 121 106)
"A sequence of random integers ranging between -255 and 255.")
(defun make-vector-from-pool (type length &optional (starting 0))
"Make a vector of the specified element type and length using the
pool data for the type and starting at the specified point in the pool."
(make-marray
type
:dimensions length
:initial-contents
(make-list-from-pool type length starting)))
(defun make-list-from-pool (type length &optional (starting 0))
"Make a list for :initial-contents of the specified element type and
length using the pool data for the type and starting at the
......
;; Macros to interface GSL functions, including definitions necessary for defmfun.
;; Liam Healy
;; Time-stamp: <2009-04-25 23:31:38EDT interface.lisp>
;; $Id$
;; Time-stamp: <2009-11-12 18:11:41EST interface.lisp>
(in-package :gsl)
......@@ -78,7 +77,7 @@
(member value (list +success+ +continue+)))
(defun success-continue (value)
"If status is +success+, return T, othewise return NIL."
"If status is +success+, return T, otherwise return NIL."
(eql value +success+))
(defun check-gsl-status (status-code context)
......
;; Number types used by GSL functions, and specification conversion
;; Liam Healy 2008-12-31 21:06:34EST types.lisp
;; Time-stamp: <2009-05-25 16:34:12EDT types.lisp>
;; Time-stamp: <2009-11-14 10:54:53EST types.lisp>
;; $Id: $
(in-package :gsl)
......@@ -207,16 +207,3 @@
"The CL type from the CFFI element type."
(unless (eq cffi-type :pointer)
(lookup-type cffi-type *cstd-cl-type-mapping*)))
(defun splice-name (base-name type keyword)
"Make a new C name for a data function from a base name."
(let ((pos (search keyword base-name)))
(when pos
(let ((insert (+ pos (length keyword)))
(gsltype (cl-gsl type)))
(concatenate 'string
(subseq base-name 0 insert)
(if (zerop (length gsltype)) "" "_")
(cl-gsl type)
(subseq base-name insert))))))
;; Evaluation of interpolation functions.
;; Liam Healy Sun Nov 4 2007 - 18:40
;; Time-stamp: <2009-01-28 19:22:15EST evaluation.lisp>
;; Time-stamp: <2009-11-11 23:18:50EST evaluation.lisp>
(in-package :gsl)
(defmfun evaluate ((interpolation interpolation) x &key xa ya acceleration)
(defmfun evaluate
((interpolation interpolation) x
&key xa ya (acceleration (make-acceleration)))
"gsl_interp_eval"
(((mpointer interpolation) :pointer) (xa :pointer) (ya :pointer) (x :double)
((mpointer acceleration) :pointer))
......@@ -24,7 +26,8 @@
xa and ya and the accelerator acceleration."))
(defmfun evaluate-derivative
((interpolation interpolation) x &key xa ya acceleration)
((interpolation interpolation) x
&key xa ya (acceleration (make-acceleration)))
"gsl_interp_eval_deriv"
(((mpointer interpolation) :pointer) (xa :pointer) (ya :pointer) (x :double)
((mpointer acceleration) :pointer))
......@@ -40,7 +43,7 @@
xa and ya and the accelerator acceleration."))
(defmfun evaluate-second-derivative
((interpolation interpolation) x &key xa ya acceleration)
((interpolation interpolation) x &key xa ya (acceleration (make-acceleration)))
"gsl_interp_eval_deriv2"
(((mpointer interpolation) :pointer) (xa :pointer) (ya :pointer) (x :double)
((mpointer acceleration) :pointer))
......@@ -56,7 +59,8 @@
data arrays xa and ya and the accelerator 'acceleration."))
(defmfun evaluate-integral
((interpolation interpolation) low high &key xa ya acceleration)
((interpolation interpolation) low high
&key xa ya (acceleration (make-acceleration)))
"gsl_interp_eval_integ"
(((mpointer interpolation) :pointer) (xa :pointer) (ya :pointer)
(low :double) (high :double) ((mpointer acceleration) :pointer))
......@@ -65,25 +69,28 @@
:c-return :double)
;;; Spline
(defmfun evaluate ((spline spline) x &key acceleration)
(defmfun evaluate ((spline spline) x &key (acceleration (make-acceleration)))
"gsl_spline_eval"
(((mpointer spline) :pointer) (x :double) ((mpointer acceleration) :pointer))
:definition :method
:c-return :double)
(defmfun evaluate-derivative ((spline spline) x &key acceleration)
(defmfun evaluate-derivative
((spline spline) x &key (acceleration (make-acceleration)))
"gsl_spline_eval_deriv"
(((mpointer spline) :pointer) (x :double) ((mpointer acceleration) :pointer))
:definition :method
:c-return :double)
(defmfun evaluate-second-derivative ((spline spline) x &key acceleration)
(defmfun evaluate-second-derivative
((spline spline) x &key (acceleration (make-acceleration)))
"gsl_spline_eval_deriv2"
(((mpointer spline) :pointer) (x :double) ((mpointer acceleration) :pointer))
:definition :method
:c-return :double)
(defmfun evaluate-integral ((spline spline) low high &key acceleration)
(defmfun evaluate-integral
((spline spline) low high &key (acceleration (make-acceleration)))
"gsl_spline_eval_integ"
(((mpointer spline) :pointer) (low :double) (high :double)
((mpointer acceleration) :pointer))
......
;; Example spline
;; Liam Healy, Sat Nov 10 2007 - 21:18
;; Time-stamp: <2009-02-16 09:46:57EST spline-example.lisp>
;; $Id$
;; Time-stamp: <2009-11-11 23:19:51EST spline-example.lisp>
(in-package :gsl)
(defun spline-example (&optional (step 0.01d0))
"The first example in Sec. 26.7 of the GSL manual."
(let* ((acc (make-acceleration))
(xarr
(let* ((xarr
(make-marray
'double-float
:initial-contents
......@@ -22,6 +20,25 @@
collect (+ i (cos (expt i 2))))))
(spline (make-spline +cubic-spline-interpolation+ xarr yarr)))
(loop for xi from (maref xarr 0) below (maref xarr 9) by step
collect (list xi (evaluate spline xi :acceleration acc)))))
collect (list xi (evaluate spline xi)))))
(save-test interpolation (spline-example 0.1d0))
(defun evaluate-integral-example (&optional (intervals 4))
"Evaluate integral of sin(x) in interval 0-pi. sin(x) is tabulated
over a 0-2pi interval and interpolated with
+periodic-cubic-spline-interpolation+"
(let* ((xarr
(loop with step = (/ (* 2.0 pi) intervals)
for i from 0 upto intervals
collect (* i step)))
(xmarr (make-marray 'double-float :initial-contents xarr))
(ymarr
(make-marray 'double-float :initial-contents (mapcar 'sin xarr)))
(spline
(make-spline +periodic-cubic-spline-interpolation+ xmarr ymarr)))
(evaluate-integral spline 0d0
#+clisp (coerce pi 'double-float) ;; pi=3.14..L0 on clisp
#-clisp pi)))
(save-test interpolation
(spline-example 0.1d0)
(evaluate-integral-example 100))
;; Tridiagonal and Bidiagonal matrices
;; Liam Healy, Thu May 4 2006 - 15:43
;; Time-stamp: <2008-12-07 18:31:38EST diagonal.lisp>
;; Time-stamp: <2009-11-14 09:22:50EST diagonal.lisp>
;; $Id$
(in-package :gsl)
......@@ -211,3 +211,40 @@
( e_0 d_1 e_1 0 )
( 0 e_1 d_2 e_2 )
( e_3 0 e_2 d_3 )")
;;;;****************************************************************************
;;;; Example
;;;;****************************************************************************
(defun solve-tridiagonal-example (&optional (n 6))
"Solution differential equation
y=1 with boundary conditions y(0)=(n-1)^2, y(n-1)=0.
The solution is the sequence: (n-1)^2, (n-2)^2, ... 9, 4, 2, 1, 0
Desicretization of y leads to a tridiagonal system of equations
(y_(i-1)-2_i+y_(i+1))/2 = 1 for 1 < i < (n - 2)
The boundary conditions are implemented as
y(0)=(n-1)^2
y(n-1)=0"
(labels ((mvec (dimension &key initial-element)
"Shorthand for initializing gsll vectors"
(if initial-element
(make-marray 'double-float :dimensions dimension :initial-element initial-element)
(make-marray 'double-float :dimensions dimension))))
(let ((x (mvec n))
(b (mvec n :initial-element 1d0))
(diag (mvec n :initial-element -1d0))
(superdiag (mvec (1- n) :initial-element 0.5d0))
(subdiag (mvec (1- n) :initial-element 0.5d0)))
(setf
;; i=0 matrix elements
(maref diag 0) 1d0
(maref superdiag 0) 0d0
(maref b 0) (coerce (expt (1- n) 2) 'double-float)
;; i=n-1 matrix elements
(maref diag (1- n)) 1d0
(maref subdiag (- n 2)) 0d0
(maref b (1- n)) 0d0)
(solve-tridiagonal diag superdiag subdiag b x))))
;; Dilogarithm
;; Liam Healy, Fri Mar 17 2006 - 18:44
;; Time-stamp: <2008-10-21 22:57:50EDT dilogarithm.lisp>
;; $Id$
;; Time-stamp: <2009-11-14 10:27:52EST dilogarithm.lisp>
(in-package :gsl)
;;; dilog merge complex and real
(export 'dilogarithm)
(defgeneric dilogarithm (x)
(:documentation ; FDL
"The dilogarithm."))
......
;; Psi (digamma) functions
;; Liam Healy, Mon May 1 2006 - 22:11
;; Time-stamp: <2008-10-25 11:27:13EDT psi.lisp>
;; Time-stamp: <2009-11-14 10:56:27EST psi.lisp>
;; $Id$
(in-package :gsl)
......@@ -36,6 +36,7 @@
;;;; Trigamma Function
;;;;****************************************************************************
(export 'psi-1)
(defgeneric psi-1 (x)
;; FDL
(:documentation "The Trigamma function."))
......
......@@ -3,7 +3,10 @@
(in-package :gsl)
(LISP-UNIT:DEFINE-TEST INTERPOLATION
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST 1.9999999566460036d0)
(MULTIPLE-VALUE-LIST (EVALUATE-INTEGRAL-EXAMPLE 100)))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST
(LIST (LIST 0.0d0 1.0d0)
(LIST 0.1d0 1.069722939658712d0)
......
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