Commit 1d0a5ab2 authored by Liam Healy's avatar Liam Healy

Merge branch 'experimental' into cffi-fsbv

Conflicts:
	init/funcallable.lisp
	linear-algebra/blas2.lisp
	linear-algebra/blas3.lisp
	solve-minimize-fit/linear-least-squares.lisp
	special-functions/coulomb.lisp
parents 36432b9b 62d44d73
;; Expand the body of a defmfun
;; Liam Healy 2009-04-13 22:07:13EDT body-expand.lisp
;; Time-stamp: <2010-12-04 09:48:48EST body-expand.lisp>
;; Time-stamp: <2010-12-16 11:46:22EST body-expand.lisp>
;;
;; Copyright 2009, 2010 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -84,20 +84,7 @@
(defun cl-convert-form (decl)
"Generate a form that calls the appropriate converter from C/GSL to CL."
(or
#+fsbv
(if (fsbv:converter-defined-p (grid:st-actual-type decl))
`((fsbv:object ,(grid:st-symbol decl) ',(grid:st-actual-type decl))))
#-fsbv nil
(case (grid:st-actual-type decl)
(sf-result
`((val ,(grid:st-symbol decl))
(err ,(grid:st-symbol decl))))
(sf-result-e10
`((val ,(grid:st-symbol decl) 'sf-result-e10)
(e10 ,(grid:st-symbol decl))
(err ,(grid:st-symbol decl) 'sf-result-e10)))
(t `((cffi:mem-aref ,(grid:st-symbol decl) ',(grid:st-actual-type decl)))))))
(list `(fsbv:object ,(grid:st-symbol decl) ',(grid:st-actual-type decl))))
(defun values-unless-singleton (forms)
(unless (listp forms) (error "Values are not a list."))
......
;; Generate a lambda that calls the user function; will be called by callback.
;; Liam Healy
;; Time-stamp: <2011-01-10 18:19:11EST funcallable.lisp>
;; Time-stamp: <2011-10-23 20:25:24EDT funcallable.lisp>
;;
;; Copyright 2009, 2010, 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -97,7 +97,7 @@
,@(affi::delinearize-index
(affi:make-affi (value-from-dimensions argspec dimension-values))
linear-index))
`(cffi:mem-aref
`(fsbv:object
,foreign-pointer-name
',(parse-callback-argspec argspec 'element-type)
,linear-index))
......
;; Macros to interface GSL functions, including definitions necessary for defmfun.
;; Liam Healy
;; Time-stamp: <2010-06-29 19:42:14EDT interface.lisp>
;; Time-stamp: <2010-12-19 15:58:10EST interface.lisp>
;;
;; Copyright 2009 Liam M. Healy
;; Copyright 2009, 2010 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
......@@ -41,10 +41,6 @@
;;;; Declare foreign objects
;;;;****************************************************************************
(defmacro scref (size &optional (index 0))
"Reference C size(s)."
`(cffi:mem-aref ,size 'sizet ,index))
(defun wfo-declare (d cbinfo)
`(,(grid:st-symbol d)
,@(if (eq (grid:st-symbol d)
......
;; Number types used by GSL functions, and specification conversion
;; Liam Healy 2008-12-31 21:06:34EST types.lisp
;; Time-stamp: <2010-06-27 18:03:22EDT types.lisp>
;; Time-stamp: <2010-12-19 11:42:56EST types.lisp>
;;
;; Copyright 2008, 2009 Liam M. Healy
;; Copyright 2008, 2009, 2010 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
......@@ -24,13 +24,15 @@
;;;; Unsigned address types size_t
;;;;****************************************************************************
(fsbv:defconvert (sizet :constructor values))
(case (cffi:foreign-type-size :long)
(8
(push :int64 *features*)
#+fsbv (fsbv:defsynonym sizet :uint64))
#+fsbv (fsbv:define-equivalent-type sizet :uint64))
(4
(push :int32 *features*)
#+fsbv (fsbv:defsynonym sizet :uint32)))
#+fsbv (fsbv:define-equivalent-type sizet :uint32)))
;;;;****************************************************************************
;;;; Type specification conversion
......
;; BLAS level 2, Matrix-vector operations
;; Liam Healy, Wed Apr 26 2006 - 21:08
;; Time-stamp: <2011-04-23 17:15:33EDT blas2.lisp>
;; Time-stamp: <2011-10-23 20:26:44EDT blas2.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy
;; Copyright 2006, 2007, 2008, 2009, 2010, 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
......@@ -17,7 +17,6 @@
;;
;; You should have received a copy of the GNU General Public License
;; along with this program. If not, see <http://www.gnu.org/licenses/>.
;; $Id$
(in-package :gsl)
......@@ -27,14 +26,12 @@
;;;; Options
;;;;****************************************************************************
#+fsbv
(fsbv:defcenum-aux cblas-transpose)
#+fsbv (fsbv:defenumeration cblas-transpose)
;;; Eventually, #-fsbv (cffi:defconvert-cenum cblas-transpose)
#+fsbv
(fsbv:defcenum-aux cblas-uplo)
#+fsbv (fsbv:defenumeration cblas-uplo)
#+fsbv
(fsbv:defcenum-aux cblas-diag)
#+fsbv (fsbv:defenumeration cblas-diag)
;;;;****************************************************************************
;;;; Functions
......
;; BLAS level 3, Matrix-matrix operations
;; Liam Healy, Wed Apr 26 2006 - 21:08
;; Time-stamp: <2011-01-12 00:49:27EST blas3.lisp>
;; Time-stamp: <2011-10-23 20:27:45EDT blas3.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy
;; Copyright 2006, 2007, 2008, 2009, 2010, 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
......@@ -26,8 +26,8 @@
;;;; Options
;;;;****************************************************************************
#+fsbv
(fsbv:defcenum-aux cblas-side)
#+fsbv (fsbv:defenumeration cblas-side)
;;; Eventually, #-fsbv (cffi:defconvert-cenum cblas-side)
;;;;****************************************************************************
;;;; Functions
......
;; Generators of random numbers.
;; Liam Healy, Sat Jul 15 2006 - 14:43
;; Time-stamp: <2010-07-16 17:11:25EDT generators.lisp>
;; Time-stamp: <2010-12-19 16:11:01EST generators.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Copyright 2006, 2007, 2008, 2009, 2010 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
......@@ -169,7 +169,7 @@
(loop for i from 0 below (length ans)
do
(setf (aref ans i)
(mem-aref (rng-state gen) :uint8 i)))
(fsbv:object (rng-state gen) :uint8 i)))
ans))
;;;;****************************************************************************
......
;; Linear least squares, or linear regression
;; Liam Healy <2008-01-21 12:41:46EST linear-least-squares.lisp>
;; Time-stamp: <2011-08-16 23:46:14EDT linear-least-squares.lisp>
;; Time-stamp: <2011-10-23 20:31:38EDT linear-least-squares.lisp>
;;
;; Copyright 2008, 2009, 2011 Liam M. Healy
;; Copyright 2008, 2009, 2010, 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
......@@ -202,7 +202,7 @@
:inputs (model weight observations)
:outputs (parameters covariance)
:switch (weight)
:return (parameters covariance (grid:dcref chisq))
:return (parameters covariance (fsbv:object chisq :double))
:export nil
:index linear-mfit
:documentation ; FDL
......@@ -247,7 +247,7 @@
:inputs (model weight observations)
:outputs (parameters covariance)
:switch (weight)
:return ((grid:dcref chisq) (scref rank))
:return ((fsbv:object chisq :double) (fsbv:object rank 'sizet))
:export nil
:index linear-mfit
:documentation ; FDL
......
;; Coulumb functions
;; Liam Healy, Sat Mar 18 2006 - 23:23
;; Time-stamp: <2011-01-10 17:59:32EST coulomb.lisp>
;; Time-stamp: <2011-10-23 20:32:16EDT coulomb.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy
;; Copyright 2006, 2007, 2008, 2009, 2010, 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
......@@ -52,9 +52,13 @@
(F sf-result) (Fp sf-result) (G sf-result) (Gp sf-result)
(exp-F (:pointer :double)) (exp-G (:pointer :double)))
:return
((val F) (val Fp) (val G) (val Gp)
(grid:dcref exp-F) (grid:dcref exp-G)
(err F) (err Fp) (err G) (err Gp))
((let ((vl (multiple-value-list (values-with-errors F Fp G Gp))))
(values-list
(append (subseq vl 0 4)
(list
(fsbv:object exp-F :double)
(fsbv:object exp-G :double))
(subseq vl 4)))))
:documentation ; FDL
"The Coulomb wave functions F_L(\eta,x),
G_{L-k}(\eta,x) and their derivatives F'_L(\eta,x), G'_{L-k}(\eta,x)
......@@ -72,7 +76,7 @@
((L-min :double) ((1- (dim0 fc-array)) :int) (eta :double) (x :double)
((grid:foreign-pointer fc-array) :pointer) (F-exponent (:pointer :double)))
:outputs (fc-array)
:return (fc-array (grid:dcref F-exponent))
:return (fc-array (fsbv:object F-exponent :double))
:documentation ; FDL
"The Coulomb wave function F_L(\eta,x) for
L = Lmin ... Lmin + kmax, storing the results in fc-array.
......@@ -89,7 +93,10 @@
((grid:foreign-pointer fc-array) :pointer) ((grid:foreign-pointer gc-array) :pointer)
(F-exponent (:pointer :double)) (G-exponent (:pointer :double)))
:outputs (fc-array gc-array)
:return (fc-array gc-array (grid:dcref F-exponent) (grid:dcref G-exponent))
:return (fc-array
gc-array
(fsbv:object F-exponent :double)
(fsbv:object G-exponent :double))
:documentation ; FDL
"The functions F_L(\eta,x),
G_L(\eta,x) for L = Lmin ... Lmin + kmax storing the
......@@ -114,7 +121,7 @@
:outputs (fc-array fcp-array gc-array gcp-array)
:return
(fc-array fcp-array gc-array gcp-array
(grid:dcref F-exponent) (grid:dcref G-exponent))
(fsbv:object F-exponent :double) (fsbv:object G-exponent :double))
:documentation ; FDL
"The functions F_L(\eta,x),
G_L(\eta,x) and their derivatives F'_L(\eta,x),
......@@ -130,7 +137,7 @@
((L-min :double) ((1- (dim0 fc-array)) :int) (eta :double) (x :double)
((grid:foreign-pointer fc-array) :pointer) (F-exponent (:pointer :double)))
:outputs (fc-array)
:return (fc-array (grid:dcref F-exponent))
:return (fc-array (fsbv:object F-exponent :double))
:documentation ; FDL
"The Coulomb wave function divided by the argument
F_L(\eta, x)/x for L = Lmin ... Lmin + kmax, storing the
......
;; Dilogarithm
;; Liam Healy, Fri Mar 17 2006 - 18:44
;; Time-stamp: <2009-12-27 10:10:05EST dilogarithm.lisp>
;; Time-stamp: <2010-12-04 11:51:29EST dilogarithm.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Copyright 2006, 2007, 2008, 2009, 2010 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
......@@ -34,7 +34,7 @@
"gsl_sf_complex_dilog_e"
(((abs x) :double) ((phase x) :double) (re sf-result) (im sf-result))
:definition :method
:return ((complex (val re) (val im)) (complex (err re) (err im))))
:return ((complex-with-error re im)))
(save-test dilogarithm
(dilogarithm 1.0d0)
......
;; Gamma functions
;; Liam Healy, Thu Apr 27 2006 - 22:06
;; Time-stamp: <2010-06-27 18:03:15EDT gamma.lisp>
;; Time-stamp: <2010-12-04 11:39:55EST gamma.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Copyright 2006, 2007, 2008, 2009, 2010 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
......@@ -52,7 +52,9 @@
(defmfun log-gamma-sign (x)
"gsl_sf_lngamma_sgn_e"
((x :double) (ret sf-result) (sign (:pointer :double)))
:return ((val ret) (grid:dcref sign) (err ret))
:return ((multiple-value-bind (val err)
(fsbv:object ret 'sf-result)
(values val (fsbv:object sign :double) err )))
:documentation ; FDL
"Compute the sign of the gamma function and the logarithm of
its magnitude, subject to x not being a negative integer. The
......@@ -89,7 +91,7 @@
will result in a :ELOSS error when it occurs. The absolute
value part (lnr), however, never suffers from loss of precision."
:return
((val lnr) (val arg) (err lnr) (err arg)))
((values-with-errors lnr arg)))
(defmfun taylor-coefficient (n x)
"gsl_sf_taylorcoeff_e" ((n :int) (x :double) (ret sf-result))
......@@ -151,7 +153,9 @@
The computed parameters are result =
log(|(a)_x|) and sgn = sgn((a)_x) where (a)_x :=
Gamma(a + x)/Gamma(a), subject to a, a+x not being negative integers."
:return ((val ret) (grid:dcref sign) (err ret)))
:return ((multiple-value-bind (val err)
(fsbv:object ret 'sf-result)
(values val (fsbv:object sign :double) err ))))
(defmfun relative-pochammer (a x)
"gsl_sf_pochrel_e" ((a :double) (x :double) (ret sf-result))
......
;; Logarithm
;; Liam Healy, Sun Apr 30 2006 - 22:08
;; Time-stamp: <2009-12-27 10:10:01EST logarithm.lisp>
;; Time-stamp: <2010-12-04 11:41:14EST logarithm.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -35,8 +35,7 @@
(((realpart x) :double) ((imagpart x) :double)
(re-ret sf-result) (im-ret sf-result))
:definition :method
:return
((complex (val re-ret) (val im-ret)) (complex (err re-ret) (err im-ret)))
:return ((complex-with-error re-ret im-ret))
:documentation ; FDL
"Results are returned as lnr, theta such that
exp(lnr + i \theta) = z_r + i z_i, where theta lies in the range [-\pi,\pi].")
......
;; Structures returned by special functions.
;; Liam Healy, Mon Jan 1 2007 - 11:35
;; Time-stamp: <2010-06-29 22:15:22EDT return-structures.lisp>
;; Time-stamp: <2010-12-19 16:09:12EST return-structures.lisp>
;;
;; Copyright 2007, 2008, 2009 Liam M. Healy
;; Copyright 2007, 2008, 2009, 2010 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
......@@ -24,6 +24,11 @@
;;;; Result from special functions
;;;;****************************************************************************
#|
;;; Remove #'val, #'err, #'e10 in favor of using #'fsbv:object to
;;; convert from foreign to CL.
(defun val (sf-result &optional (type 'sf-result))
(cffi:foreign-slot-value sf-result type 'val))
......@@ -32,6 +37,43 @@
(defun e10 (sf-result)
(cffi:foreign-slot-value sf-result 'sf-result-e10 'e10))
|#
;;; Define instead with fsbv:defconvert. Eventually this should get
;;; into CFFI and automatically from cffi:defcstruct; and also picked
;;; up by cffi-grovel. For the time being, FSBV is required.
#+fsbv
(fsbv:defconvert (sf-result :constructor values)
(val :double)
(err :double))
(defun values-e10 (val err e10)
(values val e10 err))
#+fsbv
(fsbv:defconvert (sf-result-e10 :constructor values-e10)
(val :double)
(err :double)
(e10 :int))
(defun complex-with-error (real-sfr imaginary-sfr)
"Return two complex numbers, the value and the error."
(multiple-value-bind (re-val re-err)
(fsbv:object real-sfr 'sf-result)
(multiple-value-bind (im-val im-err)
(fsbv:object imaginary-sfr 'sf-result)
(values
(complex re-val im-val)
(complex re-err im-err)))))
(defun values-with-errors (&rest values-sfr)
"Return numbers as values and errors."
(loop for vs in values-sfr
for (val err) = (multiple-value-list (fsbv:object vs 'sf-result))
collect val into values
collect err into errors
finally (return (values-list (append values errors)))))
;;;;****************************************************************************
;;;; Array construction
......
;; Trigonometry
;; Liam Healy, Thu May 4 2006 - 22:58
;; Time-stamp: <2009-12-27 10:09:58EST trigonometry.lisp>
;; Time-stamp: <2010-12-04 11:32:57EST trigonometry.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -40,8 +40,7 @@
(((realpart x) :double) ((imagpart x) :double)
(re-ret sf-result) (im-ret sf-result))
:definition :method
:return ((complex (val re-ret) (val im-ret))
(complex (err re-ret) (err im-ret))))
:return ((complex-with-error re-ret im-ret)))
(defmfun gsl-cos ((x float))
"gsl_sf_cos_e" ((x :double) (ret sf-result))
......@@ -53,8 +52,7 @@
(((realpart x) :double) ((imagpart x) :double)
(re-ret sf-result) (im-ret sf-result))
:definition :method
:return ((complex (val re-ret) (val im-ret))
(complex (err re-ret) (err im-ret))))
:return ((complex-with-error re-ret im-ret)))
(defmfun hypotenuse (x y)
"gsl_sf_hypot_e" ((x :double) (y :double) (ret sf-result))
......@@ -74,8 +72,7 @@
"This function computes the logarithm of the complex sine,
\log(\sin(z_r + i z_i)) storing the real and imaginary parts in
szr, szi."
:return ((complex (val re-ret) (val im-ret))
(complex (err re-ret) (err im-ret))))
:return ((complex-with-error re-ret im-ret)))
;;;;****************************************************************************
;;;; Hyperbolic Trigonometric Functions
......@@ -100,7 +97,7 @@
(defmfun polar-to-rectangular (r theta)
"gsl_sf_polar_to_rect"
((r :double) (theta :double) (x sf-result) (y sf-result))
:return ((val x) (val y) (err x) (err y))
:return ((values-with-errors x y))
:documentation ; FDL
"Convert the polar coordinates (r, theta) to
rectilinear coordinates (x, y), x = r\cos(\theta), y = r\sin(\theta).")
......@@ -108,7 +105,7 @@
(defmfun rectangular-to-polar (x y)
"gsl_sf_rect_to_polar"
((x :double) (y :double) (r sf-result) (theta sf-result))
:return ((val r) (val theta) (err r) (err theta))
:return ((values-with-errors r theta))
:documentation ; FDL
"Convert the rectilinear coordinates (x, y) to
polar coordinates (r, theta), such that x =
......
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