From 1ab6ee5ad43a92228c7bc462c488c4ee6a153d41 Mon Sep 17 00:00:00 2001 From: liam <liam@a3d8a0fb-c1db-0310-ace7-a616afeb9e30> Date: Sat, 29 Apr 2006 02:52:49 +0000 Subject: [PATCH] Allow arbitrary CL function arglist in defun-gsl with :function argument, like :method argument. The arguments to the GSL function can be mapped arbitrarily, and different pieces like the size of a vector (with new function #'dim0) or the pointer #'gsl-array can be put in the correct place in the arglist explicitly. Thus I can simplify #'splice-arguments and eliminate the hackish automatic insertion of the vector dimension after (or maybe before) when declared (:double *). Ported all half-dozen or so uses to the new form. Add :after cl-invalidate where needed. git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3054 a3d8a0fb-c1db-0310-ace7-a616afeb9e30 --- gsll.asd | 3 +- init/interface.lisp | 58 +++++++++-------------------- linear-algebra/qr.lisp | 56 ++++++++++++++++++++-------- polynomial.lisp | 12 ++++-- sorting.lisp | 12 ++++-- special-functions/bessel.lisp | 14 ++++--- special-functions/gegenbauer.lisp | 62 +++++++++++++++++++++++++++++++ 7 files changed, 149 insertions(+), 68 deletions(-) create mode 100644 special-functions/gegenbauer.lisp diff --git a/gsll.asd b/gsll.asd index 2143b84b..24c63f74 100644 --- a/gsll.asd +++ b/gsll.asd @@ -3,7 +3,7 @@ ; description: Definition of GSLL system ; date: ; author: Liam Healy -; modified: Fri Apr 28 2006 - 17:24 +; modified: Fri Apr 28 2006 - 22:35 ;******************************************************** ;;; $Id: $ @@ -66,4 +66,5 @@ (:file "exponential-integrals") (:file "fermi-dirac") (:file "gamma") + (:file "gegenbauer") )))) diff --git a/init/interface.lisp b/init/interface.lisp index d89c0dbb..9363e9c9 100644 --- a/init/interface.lisp +++ b/init/interface.lisp @@ -3,7 +3,7 @@ ; description: Macros to interface GSL functions. ; date: Mon Mar 6 2006 - 22:35 ; author: Liam M. Healy -; modified: Fri Apr 21 2006 - 09:26 +; modified: Fri Apr 28 2006 - 22:31 ;******************************************************** (in-package :gsl) @@ -203,46 +203,23 @@ and a scaling exponent e10, such that the value is val*10^e10." ;;; All arrays and vectors are pass with a gsl-data object that is made ;;; outside the function, and pieces are spliced in the function call. -;;; Both raw arrays and GSL structures are be passed as gsl-data objects, -;;; but their declarations are be different: -;;; Raw array (xyz (:double *)) -;;; GSL struct: (xyz gsl-vector-c) -;;; They are expanded in the foreign-funcall argument list differently: -;;; Raw array to two arguments (:POINTER (gsl-array xyz) :SIZE (storage-size xyz)) -;;; GSL struct to one: (:pointer (pointer xyz)) - -;;; (splice-arguments '((argraw (:double *)) (argstruct gsl-vector-c) (x :double))) +;;; If a raw array is needed, use the :function argument to defun-gsl, +;;; then use the appropriate function(s) e.g. gsl-array, dim0 +;;; on that array to map to the GSL arguments. (defun splice-arguments (arguments &optional mode) "Convert the argument declarations to a list of declarations appropriate for foreign-funcall. If mode is T, a mode argument will be added to the end, if it is an integer, it will be put at that position." - (let (sizes) ; save the dimension arg if not '* - (flet ((splicearg (spec) - ;; No accomodation for matrices (two indices) yet - (if (rst-arrayp spec) - (let* ((value `(first (storage-size ,(rst-symbol spec)))) - (use-form - (if (eq (rst-dim spec) '*) - value - (progn - (push `(,(rst-dim spec) ,value) sizes) - (rst-dim spec))))) - (if (rst-pointer-last-p spec) - `(:size ,use-form - :pointer (gsl-array ,(rst-symbol spec))) - `(:pointer (gsl-array ,(rst-symbol spec)) - :size ,use-form))) - `(,(rst-type spec) - ,(wrap-arg spec))))) - (values - (mapcan #'splicearg - (if mode - (let ((mode (if (integerp mode) mode (length arguments)))) - (append (subseq arguments 0 mode) - '((mode sf-mode)) - (subseq arguments mode))) - arguments)) - sizes)))) + (flet ((splicearg (spec) + `(,(rst-type spec) + ,(wrap-arg spec)))) + (mapcan #'splicearg + (if mode + (let ((mode (if (integerp mode) mode (length arguments)))) + (append (subseq arguments 0 mode) + '((mode sf-mode)) + (subseq arguments mode))) + arguments)))) ;;;;**************************************************************************** ;;;; Checking results @@ -283,8 +260,8 @@ and a scaling exponent e10, such that the value is val*10^e10." ;;; Warning isn't quite right for lambdas. (defmacro defun-gsl (cl-name arguments gsl-name - &key documentation return mode (c-return-value :error-code) - return-input check-null-pointers method after) + &key documentation return mode (c-return-value :error-code) + return-input check-null-pointers function method after) "Define a CL function that provides an interface to a GSL function. If cl-name is :lambda, make a lambda. Arguments: arguments: a list of input arguments (symbol type) to the GSL function @@ -302,9 +279,10 @@ and a scaling exponent e10, such that the value is val*10^e10." method Make output a defmethod with the value as the arglist; 'arguments should then include explicit mapping of all arguments to GSL form. + function Arguments for CL function (like :method, but make a function) after Functions to call after the GSL function has been called; result is discarded." - (let ((clargs (or method (mapcar #'rst-symbol arguments))) + (let ((clargs (or function method (mapcar #'rst-symbol arguments))) (return-symb-type (unless (or (eq c-return-value :return) return-input) (return-symbol-type return)))) diff --git a/linear-algebra/qr.lisp b/linear-algebra/qr.lisp index 4284cb65..ea78e04e 100644 --- a/linear-algebra/qr.lisp +++ b/linear-algebra/qr.lisp @@ -3,7 +3,7 @@ ; description: QR decomposition ; date: Fri Apr 28 2006 - 16:53 ; author: Liam Healy -; modified: Fri Apr 28 2006 - 17:23 +; modified: Fri Apr 28 2006 - 21:03 ;******************************************************** ;;; $Id: $ @@ -38,22 +38,28 @@ as used by @sc{lapack}. The algorithm used to perform the decomposition is Householder QR (Golub - & Van Loan, @cite{Matrix Computations}, Algorithm 5.2.1).") + & Van Loan, @cite{Matrix Computations}, Algorithm 5.2.1)." + :after ((cl-invalidate A)) + :return-input (A tau)) (defun-gsl QR-solve ((QR gsl-matrix-c) (tau gsl-vector-c) (b gsl-vector-c) (x gsl-vector-c)) "gsl_linalg_QR_solve" :documentation "Solve the square system @math{A x = b} using the @math{QR} decomposition of @math{A} into (@var{QR}, @var{tau}) given by - @code{gsl_linalg_QR_decomp}. The least-squares solution for rectangular systems can - be found using QR-lssolve.") + QR-decomp. The least-squares solution for rectangular systems can + be found using QR-lssolve." + :after ((cl-invalidate x)) + :return-input (x)) (defun-gsl QR-svx ((QR gsl-matrix-c) (tau gsl-vector-c) (x gsl-vector-c)) "gsl_linalg_QR_svx" :documentation "Solves the square system @math{A x = b} in-place using the @math{QR} decomposition of @math{A} into (@var{QR},@var{tau}) given by QR-decomp. On input @var{x} should contain the - right-hand side @math{b}, which is replaced by the solution on output.") + right-hand side @math{b}, which is replaced by the solution on output." + :after ((cl-invalidate x)) + :return-input (x)) (defun-gsl QR-lssolve ((QR gsl-matrix-c) (tau gsl-vector-c) (b gsl-vector-c) (x gsl-vector-c) @@ -65,7 +71,9 @@ residual, @math{||Ax - b||}.The routine uses the @math{QR} decomposition of @math{A} into (@var{QR}, @var{tau}) given by @code{gsl_linalg_QR_decomp}. The solution is returned in @var{x}. The - residual is computed as a by-product and stored in @var{residual}.") + residual is computed as a by-product and stored in @var{residual}." + :after ((cl-invalidate x)) + :return-input (x)) (defun-gsl QR-QTvec ((QR gsl-matrix-c) (tau gsl-vector-c) (v gsl-vector-c)) "gsl_linalg_QR_QTvec" @@ -73,7 +81,9 @@ (@var{QR},@var{tau}) to the vector @var{v}, storing the result @math{Q^T v} in @var{v}. The matrix multiplication is carried out directly using the encoding of the Householder vectors without needing to form the full - matrix @math{Q^T}.") + matrix @math{Q^T}." + :after ((cl-invalidate v)) + :return-input (v)) (defun-gsl QR-Qvec ((QR gsl-matrix-c) (tau gsl-vector-c) (v gsl-vector-c)) "gsl_linalg_QR_Qvec" @@ -81,13 +91,17 @@ (@var{QR},@var{tau}) to the vector @var{v}, storing the result @math{Q v} in @var{v}. The matrix multiplication is carried out directly using the encoding of the Householder vectors without needing to form the full - matrix @math{Q}.") + matrix @math{Q}." + :after ((cl-invalidate v)) + :return-input (v)) (defun-gsl QR-Rsolve ((QR gsl-matrix-c) (b gsl-vector-c) (x gsl-vector-c)) "gsl_linalg_QR_Rsolve" :documentation "Solve the triangular system @math{R x = b} for @var{x}. It may be useful if the product @math{b' = Q^T b} has already - been computed using QR-QTvec}.") + been computed using QR-QTvec}." + :after ((cl-invalidate x)) + :return-input (x)) (defun-gsl QR-Rsvx ((QR gsl-matrix-c) (x gsl-vector-c)) "gsl_linalg_QR_Rsvx" @@ -95,21 +109,27 @@ in-place. On input @var{x} should contain the right-hand side @math{b} and is replaced by the solution on output. This function may be useful if the product @math{b' = Q^T b} has already been computed using - QR-QTvec}.") + QR-QTvec}." + :after ((cl-invalidate x)) + :return-input (x)) (defun-gsl QR-unpack ((QR gsl-matrix-c) (tau gsl-vector-c) (Q gsl-matrix-c) (R gsl-matrix-c)) "gsl_linalg_QR_unpack" :documentation "Unpack the encoded @math{QR} decomposition (@var{QR},@var{tau}) into the matrices @var{Q} and @var{R}, where - @var{Q} is @math{M}-by-@math{M} and @var{R} is @math{M}-by-@math{N}.") + @var{Q} is @math{M}-by-@math{M} and @var{R} is @math{M}-by-@math{N}." + :after ((cl-invalidate Q R)) + :return-input (Q R)) (defun-gsl QR-QRsolve ((Q gsl-matrix-c) (R gsl-matrix-c) (b gsl-vector-c) (x gsl-vector-c)) "gsl_linalg_QR_QRsolve" :documentation "Solves the system @math{R x = Q^T b} for @var{x}. It can be used when the @math{QR} decomposition of a matrix is available in - unpacked form as (@var{Q}, @var{R}).") + unpacked form as (@var{Q}, @var{R})." + :after ((cl-invalidate x)) + :return-input (x)) (defun-gsl QR-update ((Q gsl-matrix-c) (R gsl-matrix-c) (w gsl-vector-c) (v gsl-vector-c)) @@ -118,15 +138,21 @@ decomposition (@var{Q}, @var{R}). The update is given by @math{Q'R' = Q R + w v^T} where the output matrices @math{Q'} and @math{R'} are also orthogonal and right triangular. Note that @var{w} is destroyed by the - update.") + update." + :after ((cl-invalidate w Q R)) + :return-input (Q R)) (defun-gsl R-solve ((R gsl-matrix-c) (b gsl-vector-c) (x gsl-vector-c)) "gsl_linalg_R_solve" :documentation "Solves the triangular system @math{R x = b} for the - @math{N}-by-@math{N} matrix @var{R}.") + @math{N}-by-@math{N} matrix @var{R}." + :after ((cl-invalidate x)) + :return-input (x)) (defun-gsl R-svx ((R gsl-matrix-c) (x gsl-vector-c)) "gsl_linalg_R_svx" :documentation "Solve the triangular system @math{R x = b} in-place. On input @var{x} should contain the right-hand side @math{b}, which is - replaced by the solution on output.") + replaced by the solution on output." + :after ((cl-invalidate x)) + :return-input (x)) diff --git a/polynomial.lisp b/polynomial.lisp index effd66bb..026f92fb 100644 --- a/polynomial.lisp +++ b/polynomial.lisp @@ -3,7 +3,7 @@ ; description: Polynomials ; date: Tue Mar 21 2006 - 18:33 ; author: Liam M. Healy -; modified: Sat Apr 22 2006 - 15:29 +; modified: Fri Apr 28 2006 - 22:27 ;******************************************************** ;;; $Id: $ @@ -21,8 +21,10 @@ ;;; (setf (data vec) #(1.0d0 2.0d0 3.0d0)) ;;; (polynomial-eval vec -1.0d0) ;;; 2.0d0 -(defun-gsl polynomial-eval ((coefficients (:double *)) (x :double)) +(defun-gsl polynomial-eval + (((gsl-array coefficients) :pointer) ((dim0 coefficients) :size) (x :double)) "gsl_poly_eval" + :function (coefficients x) :documentation "Evaluate the polyonomial with coefficients at the point x." :return (:double) @@ -211,9 +213,11 @@ ,workspace)))) (defun-gsl polynomial-solve-ws - ((coefficients (:double n)) (workspace poly-complex-workspace)) + (((gsl-array coefficients) :pointer) ((dim0 coefficients) :size) + (workspace poly-complex-workspace)) "gsl_poly_complex_solve" - :return ((gsl-complex (1- n)))) + :function (coefficients workspace) + :return ((gsl-complex (1- (dim0 coefficients))))) #+future (export '(polynomial-solve)) diff --git a/sorting.lisp b/sorting.lisp index 7bce331e..bab9243a 100644 --- a/sorting.lisp +++ b/sorting.lisp @@ -3,7 +3,7 @@ ; description: Sorting ; date: Fri Apr 14 2006 - 20:20 ; author: Liam M. Healy -; modified: Wed Apr 19 2006 - 17:43 +; modified: Fri Apr 28 2006 - 22:21 ;******************************************************** ;;; $Id: $ @@ -74,20 +74,26 @@ place. The first element of @var{p} gives the index of the least element in @var{v}, and the last element of @var{p} gives the index of the greatest element in @var{v}. The vector @var{v} is not changed.") -(defun-gsl sort-vector-smallest ((dest (:double *)) (v gsl-vector-c)) +(defun-gsl sort-vector-smallest + (((gsl-array dest) :pointer) ((dim0 dest) :size) (v gsl-vector-c)) "gsl_sort_vector_smallest" + :function (dest v) :documentation "Find the smallest elements of the vector @var{v} and put them into dest, which must be shorter than v." :c-return-value :void + :after ((cl-invalidate dest)) :return-input (dest)) -(defun-gsl sort-vector-largest ((dest (:double *)) (v gsl-vector-c)) +(defun-gsl sort-vector-largest + (((gsl-array dest) :pointer) ((dim0 dest) :size) (v gsl-vector-c)) "gsl_sort_vector_largest" + :function (dest v) :documentation "Find the largest elements of the vector @var{v} and put them into dest, which must be shorter than v." :c-return-value :void + :after ((cl-invalidate dest)) :return-input (dest)) ;;;;**************************************************************************** diff --git a/special-functions/bessel.lisp b/special-functions/bessel.lisp index 95da6ec8..742af2fb 100644 --- a/special-functions/bessel.lisp +++ b/special-functions/bessel.lisp @@ -3,7 +3,7 @@ ; description: Bessel functions ; date: Fri Mar 17 2006 - 18:42 ; author: Liam M. Healy -; modified: Wed Apr 19 2006 - 17:06 +; modified: Fri Apr 28 2006 - 22:16 ;******************************************************** (in-package :gsl) @@ -321,20 +321,24 @@ "The regular cylindrical Bessel function of fractional order @math{\nu}, @math{J_\nu(x)}." :return (sf-result)) -;;; (double @var{nu}, gsl_mode_t @var{mode}, size_t @var{size}, double @var{v}[]) -;;; The mode argument comes before v, need to make that an option on -;;; the :mode argument to defun-gsl. Need to return v as well. -(defun-gsl bessel-sequence-Jnu ((nu :double) (v (:double * t))) +(defun-gsl bessel-sequence-Jnu + ((nu :double) ((dim0 v) :size) ((gsl-array v) :pointer)) "gsl_sf_bessel_sequence_Jnu_e" + :function (nu v) :documentation "The regular cylindrical Bessel function of fractional order @math{\nu}, @math{J_\nu(x)}, evaluated at a series of @math{x} values. The array @var{v} of length @var{size} contains the @math{x} values. They are assumed to be strictly ordered and positive. The array is over-written with the values of @math{J_\nu(x_i)}." + :after ((cl-invalidate v)) :return-input (v) :mode 1) +;;; (defparameter vec (make-data 'vector nil 3)) +;;; (setf (data vec) #(1.0d0 2.0d0 3.0d0)) +;;; (bessel-sequence-Jnu 0.5d0 vec) + ;;;;**************************************************************************** ;;;; Irregular Bessel Function - Fractional Order ;;;;**************************************************************************** diff --git a/special-functions/gegenbauer.lisp b/special-functions/gegenbauer.lisp new file mode 100644 index 00000000..7344e600 --- /dev/null +++ b/special-functions/gegenbauer.lisp @@ -0,0 +1,62 @@ +;******************************************************** +; file: gegenbauer.lisp +; description: Gegenbauer polynomials +; date: Fri Apr 28 2006 - 20:40 +; author: Liam M. Healy +; modified: Fri Apr 28 2006 - 22:44 +;******************************************************** +;;; $Id: $ + +(in-package :gsl) + +(defun-gsl gegenbauer-1 ((lambda :double) (x :double)) + "gsl_sf_gegenpoly_1_e" + :return (sf-result) + :documentation "The Gegenbauer polynomial @math{C^@{(\lambda)@}_1(x)}.") + +(defun-gsl gegenbauer-2 ((lambda :double) (x :double)) + "gsl_sf_gegenpoly_2_e" + :return (sf-result) + :documentation "The Gegenbauer polynomial @math{C^@{(\lambda)@}_2(x)}.") + +(defun-gsl gegenbauer-3 ((lambda :double) (x :double)) + "gsl_sf_gegenpoly_3_e" + :return (sf-result) + :documentation "The Gegenbauer polynomial @math{C^@{(\lambda)@}_3(x)}.") + +(defun-gsl gegenbauer ((n :int) (lambda :double) (x :double)) + "gsl_sf_gegenpoly_n_e" + :return (sf-result) + :documentation "The Gegenbauer polynomial + @math{C^@{(\lambda)@}_n(x)} for a specific value of @var{n}, + @var{lambda}, @var{x} subject to @math{\lambda > -1/2}, + @math{n >= 0}.") + +(defun-gsl gegenbauer-array + (((dim0 result) :int) + (lambda :double) (x :double) ((gsl-array result) :pointer)) + "gsl_sf_gegenpoly_array" + :function (lambda x result) + :documentation "Compute an array of Gegenbauer polynomials + @math{C^@{(\lambda)@}_n(x)} for @math{n = 0, 1, 2, \dots, nmax}, subject + to @math{\lambda > -1/2}, @math{nmax >= 0}." + :after ((cl-invalidate result)) + :return-input (result)) + +;;; (defparameter vec (make-data 'vector nil 3)) +;;; (gegenbauer-array 1.0d0 3.0d0 vec) + +;;;;**************************************************************************** +;;;; Examples and unit test +;;;;**************************************************************************** + +(lisp-unit:define-test gegenbauer + (lisp-unit:assert-first-fp-equal + "0.600000000000d+01" + (gegenbauer-1 1.0d0 3.0d0)) + (lisp-unit:assert-first-fp-equal + "0.350000000000d+02" + (gegenbauer-2 1.0d0 3.0d0)) + (lisp-unit:assert-first-fp-equal + "0.204000000000d+03" + (gegenbauer-3 1.0d0 3.0d0))) -- GitLab