Skip to content
Snippets Groups Projects
Commit 826fa44a authored by liam's avatar liam
Browse files

Expand defun-sf to work with arrays passed as input, and make use of

this feature in polynomial.lisp.


git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3000 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent 515a044c
No related branches found
No related tags found
No related merge requests found
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
; description: Load GSL ; description: Load GSL
; date: Sat Mar 4 2006 - 18:53 ; date: Sat Mar 4 2006 - 18:53
; author: Liam M. Healy ; author: Liam M. Healy
; modified: Fri Mar 24 2006 - 15:30 ; modified: Fri Mar 24 2006 - 20:04
;******************************************************** ;********************************************************
(defpackage gsll (defpackage gsll
...@@ -22,6 +22,5 @@ ...@@ -22,6 +22,5 @@
(cffi:use-foreign-library libgsl) (cffi:use-foreign-library libgsl)
;;; If cffi-unix is unavailable, uncomment the appropriate line: ;;; If cffi-unix is unavailable, uncomment the following line:
;;; (cffi:defctype :size :unsigned-long) ; for 64 bit computers ;;; (cffi:defctype :size :unsigned-long)
;;; (cffi:defctype :size :uint) ; for 32 bit computers
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
; description: Macros to interface GSL functions. ; description: Macros to interface GSL functions.
; date: Mon Mar 6 2006 - 22:35 ; date: Mon Mar 6 2006 - 22:35
; author: Liam M. Healy ; author: Liam M. Healy
; modified: Fri Mar 24 2006 - 17:43 ; modified: Sat Mar 25 2006 - 21:50
;******************************************************** ;********************************************************
(in-package :gsl) (in-package :gsl)
...@@ -137,12 +137,65 @@ and a scaling exponent e10, such that the value is val*10^e10." ...@@ -137,12 +137,65 @@ and a scaling exponent e10, such that the value is val*10^e10."
(gsl-complex `((complex-to-cl ,(rst-symbol decl)))) (gsl-complex `((complex-to-cl ,(rst-symbol decl))))
(:double `((double-to-cl ,(rst-symbol decl))))))) (:double `((double-to-cl ,(rst-symbol decl)))))))
(defun input-argument (spec) ;;;;****************************************************************************
"Create a CFFI input argument specification." ;;;; Wrapping for arrays
;;;;****************************************************************************
(defun wrap-form (form wrappers)
(if wrappers
(wrap-form
`(,@(first wrappers) ,form)
(rest wrappers))
form))
(defun wfa-wrapper (spec)
"Returns two values: the foreign-funcall argument declaration
and what to wrap around the form, if anything."
(if (rst-arrayp spec) (if (rst-arrayp spec)
;; currently only works for vectors (1dim array) (let ((arr (gensym "ARR"))
(list :pointer (rst-symbol spec) :size (rst-dim spec)) (len (if (eq (rst-dim spec) '*)
(list (rst-type spec) (rst-symbol spec)))) (gensym "LEN")
(rst-dim spec))))
(values
`(:pointer ,arr :size ,len)
`((cffi::with-foreign-array
(,arr ,(rst-symbol spec) ,(rst-eltype spec) (list ,len)))
(let ((,len (length ,(rst-symbol spec))))))))
(values `(,(rst-type spec) ,(rst-symbol spec)) nil)))
(defun wfa-wrappers (specs)
(if specs
(multiple-value-bind (decl wrap)
(wfa-wrapper (first specs))
(multiple-value-bind (decls wraps)
(wfa-wrappers (rest specs))
(values (append decl decls) (append wrap wraps))))
(values nil nil)))
#+example
(multiple-value-bind (decl wrap)
(wfa-wrappers
'((foo :double) (coefficients (:double *)) (n :int) (boo (:double *))))
(apply #'wrap-form `(the-body ,@decl) wrap))
#+example
(LET ((#:LEN4183 (LENGTH BOO)))
(CFFI::WITH-FOREIGN-ARRAY (#:ARR4182 BOO :DOUBLE (LIST #:LEN4183))
(LET ((#:LEN4181 (LENGTH COEFFICIENTS)))
(CFFI::WITH-FOREIGN-ARRAY
(#:ARR4180 COEFFICIENTS :DOUBLE (LIST #:LEN4181))
(THE-BODY :DOUBLE
FOO
:POINTER
#:ARR4180
:SIZE
#:LEN4181
:INT
N
:POINTER
#:ARR4182
:SIZE
#:LEN4183)))))
;;;;**************************************************************************** ;;;;****************************************************************************
;;;; Macro defun-sf ;;;; Macro defun-sf
...@@ -159,14 +212,7 @@ and a scaling exponent e10, such that the value is val*10^e10." ...@@ -159,14 +212,7 @@ and a scaling exponent e10, such that the value is val*10^e10."
if it is not :SUCCESS." if it is not :SUCCESS."
(unless (eql :success (cffi:foreign-enum-keyword 'gsl-errorno status-code)) (unless (eql :success (cffi:foreign-enum-keyword 'gsl-errorno status-code))
(warn 'gsl-warning (warn 'gsl-warning
:gsl-errno status :gsl-context context))) :gsl-errno status-code :gsl-context context)))
#+development
(defun array-input-with-dim (cl-array)
cffi:foreign-array-alloc
cffi:lisp-array-to-foreign
cffi:foreign-array-free
)
;;; Warning isn't quite right for lambdas. ;;; Warning isn't quite right for lambdas.
;;; New name? ;;; New name?
...@@ -180,42 +226,55 @@ and a scaling exponent e10, such that the value is val*10^e10." ...@@ -180,42 +226,55 @@ and a scaling exponent e10, such that the value is val*10^e10."
documentation: a string documentation: a string
return: a list of return types return: a list of return types
mode: T or NIL, depending on whether gsl_mode is an argument mode: T or NIL, depending on whether gsl_mode is an argument
c-return-value: The C function returns an :error-code or :number-of-answers." c-return-value: The C function returns an :error-code, :number-of-answers,
(let ((args (mapcar #'first arguments)) a value to :return from the CL function, or :void."
(let ((clargs (mapcar #'rst-symbol arguments))
(return-symb-type (return-symb-type
(mapcar (lambda (typ) (unless (eq c-return-value :return)
(list (gensym "RET") (mapcar (lambda (typ)
typ)) (list (gensym "RET")
return))) typ))
`(,@(if (eq cl-name :lambda) return))))
'(lambda) (multiple-value-bind (input-declarations wrap-arrays)
`(defunx-map ,cl-name ,gsl-name)) (wfa-wrappers arguments)
,(if mode `(,@(if (eq cl-name :lambda)
`(,@args &optional (mode :double-prec)) '(lambda)
`(,@args)) `(defunx-map ,cl-name ,gsl-name))
,@(when documentation (list documentation)) ,(if mode
(cffi:with-foreign-objects `(,@clargs &optional (mode :double-prec))
,(mapcar #'wfo-declare return-symb-type) `(,@clargs))
(let ((status ,@(when documentation (list documentation))
(cffi:foreign-funcall ,(wrap-form ; with-foreign-array
,gsl-name (wrap-form ; with-foreign-object
,@(mapcan #'input-argument arguments) `(let ((creturn
,@(when mode '(sf-mode mode)) (cffi:foreign-funcall
,@(mapcan (lambda (r) `(:pointer ,(rst-symbol r))) return-symb-type) ,gsl-name
:int))) ,@input-declarations
,@(if (eq c-return-value :error-code) ,@(when mode '(sf-mode mode))
`((check-gsl-status status `(,',cl-name ,,@args)))) ,@(mapcan (lambda (r) `(:pointer ,(rst-symbol r)))
(values return-symb-type)
,@(case c-return-value ,(if (eq c-return-value :return)
(:number-of-answers (first return)
(mapcan :int))))
(lambda (decl seq) ,@(case c-return-value
`((when (> status ,seq) ,@(pick-result decl)))) (:void '((declare (ignore creturn))))
return-symb-type (:error-code
(loop for i below (length return) collect i))) `((check-gsl-status creturn `(,',cl-name ,,@clargs)))))
(:return '(status)) (values
(:error-code ,@(case c-return-value
(mapcan #'pick-result return-symb-type))))))))) (:number-of-answers
(mapcan
(lambda (decl seq)
`((when (> creturn ,seq) ,@(pick-result decl))))
return-symb-type
(loop for i below (length return) collect i)))
(:return '(creturn))
(:error-code
(mapcan #'pick-result return-symb-type)))))
(when return-symb-type
`((cffi:with-foreign-objects
,(mapcar #'wfo-declare return-symb-type)))))
wrap-arrays)))))
;;; arguments: list like ((x :double) (y :double)) ;;; arguments: list like ((x :double) (y :double))
;;; mode: t or nil ;;; mode: t or nil
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
; description: Polynomials ; description: Polynomials
; date: Tue Mar 21 2006 - 18:33 ; date: Tue Mar 21 2006 - 18:33
; author: Liam M. Healy ; author: Liam M. Healy
; modified: Fri Mar 24 2006 - 15:09 ; modified: Sat Mar 25 2006 - 21:55
;******************************************************** ;********************************************************
;;; $Id: $ ;;; $Id: $
...@@ -18,16 +18,14 @@ ...@@ -18,16 +18,14 @@
;;;; Polynomial Evaluation ;;;; Polynomial Evaluation
;;;;**************************************************************************** ;;;;****************************************************************************
(defunx-map polynomial-eval "gsl_poly_eval" (coefficients x) ;;; (polynomial-eval #(1.0d0 2.0d0 3.0d0) -1.0d0)
(defun-sf polynomial-eval
((coefficients (:double *)) (x :double))
"gsl_poly_eval"
:documentation
"Evaluate the polyonomial with coefficients at the point x." "Evaluate the polyonomial with coefficients at the point x."
(let ((len (length coefficients))) :return (:double)
(cffi::with-foreign-array (coef coefficients :double (list len)) :c-return-value :return)
(cffi:foreign-funcall
"gsl_poly_eval"
:pointer coef
:int len
:double x
:double))))
;;;;**************************************************************************** ;;;;****************************************************************************
;;;; Divided Difference Representation of Polynomials ;;;; Divided Difference Representation of Polynomials
...@@ -172,7 +170,7 @@ ...@@ -172,7 +170,7 @@
(matrix :pointer)) (matrix :pointer))
(export '(with-poly-complex-workspace)) (export '(with-poly-complex-workspace))
(defmacro with-poly-complex-workspace (workspace &body body) (defmacro with-poly-complex-workspace ((workspace size) &body body)
"Macro to create and cleanup workspace for polynomial root solver." "Macro to create and cleanup workspace for polynomial root solver."
`(let ((,workspace `(let ((,workspace
(funcall (funcall
...@@ -180,16 +178,18 @@ ...@@ -180,16 +178,18 @@
"gsl_poly_complex_workspace_alloc" "gsl_poly_complex_workspace_alloc"
:return :return
(poly-complex-workspace) (poly-complex-workspace)
:c-return-value :return)))) :c-return-value :return)
,size)))
(unwind-protect (unwind-protect
(progn ,@body) (progn ,@body)
(funcall (funcall
(defun-sf :lambda ((,workspace poly-complex-workspace)) (defun-sf :lambda ((,workspace poly-complex-workspace))
"gsl_poly_complex_workspace_free" "gsl_poly_complex_workspace_free"
:c-return-value :void))))) :c-return-value :void)
,workspace))))
(defun-sf polynomial-solve (defun-sf polynomial-solve
((a :pointer) (n :size) (workspace poly-complex-workspace)) ((coefficients (:double n)) (workspace poly-complex-workspace))
"gsl_poly_complex_solve" "gsl_poly_complex_solve"
:documentation :documentation
"The roots of the general polynomial "The roots of the general polynomial
...@@ -201,4 +201,10 @@ highest order term must be non-zero. The function requires a workspace ...@@ -201,4 +201,10 @@ highest order term must be non-zero. The function requires a workspace
@var{w} of the appropriate size. The @math{n-1} roots are returned in @var{w} of the appropriate size. The @math{n-1} roots are returned in
the packed complex array @var{z} of length @math{2(n-1)}, alternating the packed complex array @var{z} of length @math{2(n-1)}, alternating
real and imaginary parts." real and imaginary parts."
:return ((gsl-complex n))) :return ((gsl-complex (1- n))))
#|
;;; Example from GSL manual
(with-poly-complex-workspace (ws 6)
(polynomial-solve #(-1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0) ws))
|#
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment