diff --git a/gsll.asd b/gsll.asd index da8432d67e3ca096808931716b869fbdc2c8af54..4cc7eccaf288d9e2b209dbbe576774bfe41c60cb 100644 --- a/gsll.asd +++ b/gsll.asd @@ -1,6 +1,6 @@ ;; Definition of GSLL system ;; Liam Healy -;; Time-stamp: <2009-04-04 22:32:18EDT gsll.asd> +;; Time-stamp: <2009-04-11 17:43:54EDT gsll.asd> ;; $Id$ (asdf:defsystem "gsll" @@ -9,7 +9,7 @@ :version "0" :author "Liam M. Healy" :licence "LLGPL v3, FDL" - :depends-on (cffi trivial-garbage) + :depends-on (cffi trivial-garbage fsbv) :components ((:module init :components diff --git a/init/complex-types.lisp b/init/complex-types.lisp index add28082f44db1c55e3c843cf15fa514e1b98d0d..c9b047133acdf5251578f5bd83e3c6d5ddb675b0 100644 --- a/init/complex-types.lisp +++ b/init/complex-types.lisp @@ -1,6 +1,6 @@ ;; Complex number types ;; Liam Healy 2009-01-13 21:24:05EST complex-types.lisp -;; Time-stamp: <2009-01-14 22:45:53EST complex-types.lisp> +;; Time-stamp: <2009-04-11 21:15:24EDT complex-types.lisp> ;; $Id: $ (in-package :gsl) @@ -10,31 +10,18 @@ ;;;;**************************************************************************** ;;; GSL represents complex variables as a struct of two reals and -;;; passes them by value (not a pointer) to functions. Few FFIs have -;;; the capability of handling that. However, it does seem to work -;;; (SBCL and CCL at least) that we can treat each complex argument as -;;; a succession of two double-float arguments, or two single-floats -;;; packed into a double float. Here we probe that the first is true; -;;; if so, functions will automatically be converted. -(defparameter *pass-complex-scalar-as-two-reals* - (eql 5.0d0 - (ignore-errors - (cffi:foreign-funcall "gsl_complex_abs" :double 3.0d0 :double 4.0d0 :double))) - "A complex number can be passed as two adjacent reals.") - -;;; GSL defines complex numbers in a struct, and passes the struct by -;;; value. CFFI does not support call by value for structs, so we -;;; cannot use functions that call or return complex scalars. +;;; passes them by value (not a pointer) to functions. In order to +;;; handle this we use FSBV. ;;; See /usr/include/gsl/gsl_complex.h -(cffi:defcstruct complex-float-c +(fsbv:defcstruct complex-float-c (dat :float :count 2)) -(cffi:defcstruct complex-double-c +(fsbv:defcstruct complex-double-c (dat :double :count 2)) #+long-double -(cffi:defcstruct complex-long-double-c +(fsbv:defcstruct complex-long-double-c (dat :long-double :count 2)) (defun clean-type (type) @@ -55,6 +42,26 @@ (second eltype) eltype)) +(defun component-float-type-from-value (value) + (let ((full (component-float-type (type-of value)))) + (if (listp full) (first full) full))) + +;; See also complex-to-cl in number-conversion.lisp +(defun cl-to-complex (complex-number struct) + "Set the CL complex number into the GSL struct." + (let ((struct-type + (ecase (component-float-type-from-value complex-number) + (single-float 'complex-float-c) + (double-float 'complex-double-c)))) + (setf (cffi:mem-aref + (cffi:foreign-slot-value struct struct-type 'dat) + :double 0) + (realpart complex-number) + (cffi:mem-aref + (cffi:foreign-slot-value struct struct-type 'dat) + :double 1) + (imagpart complex-number)))) + (define-condition pass-complex-by-value (error) () (:report @@ -65,28 +72,3 @@ "An error indicating that this implementation and platform are unable to pass complex numbers to the GSL libary by value.")) -(defun pack-complex-as-double (number) - "Pack a number of type (complex single-float) into a double-float, - which may be accepted as a struct of type gsl_complex_float. - This is hideously non-portable." - (integer-as-float - (dpb - (float-as-integer (imagpart number) t) - (byte 32 32) - (dpb - (float-as-integer (realpart number) t) - (byte 32 0) - 0)) - 'double-float)) - -(defun passing-complex-by-value (cfind) - "Substitution in defmfun so that complex numbers can be passed by - value to the GSL structs defined above. This is non-portable." - (if (eq (component-float-type (cffi-cl (third cfind))) 'single-float) - ;; single-float - `(:double (pack-complex-as-double ,(first cfind))) - ;; double-float - `(:double - (realpart ,(first cfind)) - :double - (imagpart ,(first cfind))))) diff --git a/init/defmfun-single.lisp b/init/defmfun-single.lisp index a4aff06153d74b1e1e1b9ac04b6ed17497ac010b..96a54c1e783bf4027e839bc51382490851b9fc93 100644 --- a/init/defmfun-single.lisp +++ b/init/defmfun-single.lisp @@ -1,6 +1,6 @@ ;; Helpers that define a single GSL function interface ;; Liam Healy 2009-01-07 22:02:20EST defmfun-single.lisp -;; Time-stamp: <2009-04-04 21:56:30EDT defmfun-single.lisp> +;; Time-stamp: <2009-04-11 21:24:36EDT defmfun-single.lisp> ;; $Id: $ (in-package :gsl) @@ -12,15 +12,27 @@ (second name))) (defun complex-scalars (cl-arguments c-arguments-types) + "Create a list of (symbol temp-variable c-struct) for each complex + variable in the argument." (loop for sd in c-arguments-types - for cl-type = (cffi-cl (st-type sd)) - append - (when (and cl-type (subtypep cl-type 'complex) - (member (st-symbol sd) - (arglist-plain-and-categories cl-arguments))) - (list (list (st-symbol sd) - (gensym (string (st-symbol sd))) - (st-type sd)))))) + for cl-type = (cffi-cl (st-type sd)) + append + (when (and cl-type (subtypep cl-type 'complex) + (member (st-symbol sd) + (arglist-plain-and-categories cl-arguments))) + (list (list (st-symbol sd) + (make-symbol (string (st-symbol sd))) + (st-type sd)))))) + +(defun substitute-symbols (substlist list) + "Substitute in list the second element of each list of substlist for + the first element." + (loop for pair in substlist + with answer = list + do + (setf answer + (subst (second pair) (first pair) answer)) + finally (return answer))) (defun stupid-code-walk-find-variables (sexp) "This will work with the simplest s-expression forms only to find @@ -150,83 +162,93 @@ (cret-name (if (listp c-return) (st-symbol c-return) (make-symbol "CRETURN"))) (complex-args (complex-scalars arglist c-arguments)) + (cargs (substitute-symbols complex-args c-arguments)) ; with complex (allocated ; Foreign objects to be allocated (remove-if (lambda (s) (member s (arglist-plain-and-categories arglist nil))) - (variables-used-in-c-arguments c-arguments))) + (variables-used-in-c-arguments cargs))) (allocated-decl (mapcar (lambda (s) - (or (find s c-arguments :key #'st-symbol) + (or (find s cargs :key #'st-symbol) ;; Catch programming errors, usually typos (error "Could not find ~a among the arguments" s))) allocated)) + ;; same as allocated-decl, but no complex struct variables + (allocated-decl-ret + (mapcar + (lambda (s) + (or (find s (set-difference c-arguments (mapcar 'first complex-args)) + :key #'st-symbol) + ;; Catch programming errors, usually typos + (error "Could not find ~a among the arguments" s))) + (set-difference allocated (mapcar 'second complex-args)))) (clret (or ; better as a symbol macro (substitute cret-name :c-return (mapcar (lambda (sym) - (let ((it (find sym allocated-decl :key 'st-symbol))) + (let ((it (find sym allocated-decl-ret :key 'st-symbol))) (if it (first (cl-convert-form it)) sym))) return)) (mappend #'cl-convert-form - (callback-remove-arg allocated-decl cbinfo 'st-symbol)) + (callback-remove-arg allocated-decl-ret cbinfo 'st-symbol)) outputs (unless (eq c-return :void) (list cret-name))))) - (if (and complex-args (not *pass-complex-scalar-as-two-reals*)) - '(error 'pass-complex-by-value) ; arglist should be declared ignore - (wrap-letlike - allocated-decl - (mapcar (lambda (d) (wfo-declare d cbinfo)) - allocated-decl) - 'cffi:with-foreign-objects - `(,@(append - (callback-symbol-set - callback-dynamic cbinfo (first callback-dynamic-variables)) - before - (when callback-object (callback-set-dynamic callback-object arglist))) - ,@(callback-set-slots - cbinfo callback-dynamic-variables callback-dynamic) - (let ((,cret-name - (cffi:foreign-funcall - ,gsl-name - ,@(mappend - (lambda (arg) - (let ((cfind ; variable is complex - (find (st-symbol arg) complex-args :key 'first))) - (if cfind ; make two successive scalars - (passing-complex-by-value cfind) - ;; otherwise use without conversion - (list (if (member (st-symbol arg) allocated) - :pointer - (st-type arg)) - (st-symbol arg))))) - c-arguments) - ,cret-type))) - ,@(case c-return - (:void `((declare (ignore ,cret-name)))) - (:error-code ; fill in arguments - `((check-gsl-status ,cret-name - ',(or (defgeneric-method-p name) name))))) - #-native - ,@(when outputs - (mapcar - (lambda (x) `(setf (cl-invalid ,x) t (c-invalid ,x) nil)) - outputs)) - ,@(when (eq cret-type :pointer) - `((check-null-pointer - ,cret-name - ,@'('memory-allocation-failure "No memory allocated.")))) - ,@after - (values - ,@(defmfun-return - c-return cret-name clret allocated - return return-supplied-p - enumeration outputs))))))))) + (wrap-letlike + allocated-decl + (mapcar (lambda (d) (wfo-declare d cbinfo)) + allocated-decl) + 'cffi:with-foreign-objects + `(,@(append + (callback-symbol-set + callback-dynamic cbinfo (first callback-dynamic-variables)) + before + (when complex-args + (loop for arg in complex-args + collect `(cl-to-complex ,(first arg) ,(second arg)))) + (when callback-object (callback-set-dynamic callback-object arglist))) + ,@(callback-set-slots + cbinfo callback-dynamic-variables callback-dynamic) + (let ((,cret-name + (,(if complex-args + 'fsbv:foreign-funcall + 'cffi:foreign-funcall) + ,gsl-name + ,@(mappend + (lambda (arg) + (list (cond + ((third (find (st-symbol arg) complex-args :key 'second))) + ((member (st-symbol arg) allocated) :pointer) + (t (st-type arg))) + (st-symbol arg))) + cargs) + ,cret-type))) + ,@(case c-return + (:void `((declare (ignore ,cret-name)))) + (:error-code ; fill in arguments + `((check-gsl-status ,cret-name + ',(or (defgeneric-method-p name) name))))) + #-native + ,@(when outputs + (mapcar + (lambda (x) `(setf (cl-invalid ,x) t (c-invalid ,x) nil)) + outputs)) + ,@(when (eq cret-type :pointer) + `((check-null-pointer + ,cret-name + ,@'('memory-allocation-failure "No memory allocated.")))) + ,@after + (values + ,@(defmfun-return + c-return cret-name clret + (set-difference allocated (mapcar 'second complex-args)) + return return-supplied-p + enumeration outputs)))))))) (defun defmfun-return (c-return cret-name clret allocated return return-supplied-p enumeration outputs)