error numbers and whitespace fixed.
[CommonLispStat.git] / src / numerics / ladata.lsp
blobb38c2080ef62b73b06bc1aaf66716d297528645a
1 ;;; -*- mode: lisp -*-
2 ;;; Copyright (c) 2005--2008, by A.J. Rossini <blindglobe@gmail.com>
3 ;;; See COPYRIGHT file for any additional restrictions (BSD license).
4 ;;; Since 1991, ANSI was finally finished. Edited for ANSI Common Lisp.
6 ;;;; ladata -- Data handling functions for linear algebra interface
7 ;;;;
8 ;;;; Copyright (c) 1991, by Luke Tierney. Permission is granted for
9 ;;;; unrestricted use.
11 (in-package :lisp-stat-linalg-data)
13 #+openmcl
14 (defctype size-t :unsigned-long)
15 #+sbcl
16 (defctype size-t :unsigned-int)
18 ;; Should we do the same with int's and long's? There is some
19 ;; evidence that it might be useful?
22 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
23 ;;;
24 ;;; Data Mode Functions
25 ;;;
26 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
28 ;;;;
29 ;;;; These constants need to be redefined if IN, RE or CX in linalg.h
30 ;;;; are changed.
31 ;;;;
33 ;;; FIXME:AJR: This is how Luke got around having appropriate
34 ;;; approaches for Linear Algebra. We want to cheat and instead use
35 ;;; CLEM or MATLISP as the underlying linear algebra package.
37 (defparameter +mode-in+ 0)
38 (defparameter +mode-re+ 1)
39 (defparameter +mode-cx+ 2)
41 (defun mode-of (x)
42 (etypecase x
43 (fixnum +mode-in+)
44 (rational +mode-re+)
45 (float +mode-re+)
46 (complex +mode-cx+)))
48 (defun la-data-mode (data)
49 (let ((data (compound-data-seq data))
50 (mode 0))
51 (cond
52 ((vectorp data)
53 (let ((n (length data)))
54 (declare (fixnum n))
55 (dotimes (i n mode)
56 (declare (fixnum i))
57 (setf mode (max mode (mode-of (aref data i)))))))
58 ((consp data) (dolist (x data mode) (setf mode (max mode (mode-of x)))))
59 (t (error "bad sequence - ~s" data)))))
62 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
63 ;;;;
64 ;;;; Internal Allocation Funcitons
65 ;;;;
66 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
68 (defvar *la-allocations* nil)
70 ;;;
71 ;;; CFFI glue for... Storage Allocation Functions
72 ;;;
74 (defun null-ptr-p (p) (cffi:null-pointer-p p))
75 (defun ptr-eq (p q) (cffi:pointer-eq p q))
77 (cffi:defcfun ("la_base_allocate" ccl-la-base-allocate)
78 :pointer (n size-t) (m size-t))
79 (defun la-base-allocate (n m)
80 (ccl-la-base-allocate n m))
82 (cffi:defcfun ("la_base_free_alloc" ccl-la-base-free-alloc)
83 :void (p :pointer))
84 (defun la-base-free (p)
85 (ccl-la-base-free-alloc p))
87 (cffi:defcfun ("la_mode_size" ccl-la-mode-size)
88 size-t (x size-t))
90 (defun la-mode-size (mode)
91 (ccl-la-mode-size mode))
93 ;;;
94 ;;; Callbacks for Internal Storage
95 ;;;
97 (cffi:defcallback lisp-la-allocate :void ((n size-t) (m size-t))
98 (ccl-store-ptr (la-allocate n m)))
99 (cffi:defcfun ("register_la_allocate" register-la-allocate)
100 :void (p :pointer))
101 (register-la-allocate (cffi:callback lisp-la-allocate))
102 (cffi:defcfun ("la_allocate" la)
103 :pointer (x size-t) (y size-t))
105 (cffi:defcallback lisp-la-free-alloc
106 :void ((p :pointer))
107 (la-free p))
109 (cffi:defcfun ("register_la_free_alloc" register-la-free-alloc)
110 :void (p :pointer))
111 (register-la-free-alloc (cffi:callback lisp-la-free-alloc))
112 (cffi:defcfun ("la_free_alloc" lf)
113 :void (p :pointer))
117 ;;; CFFI glue for... Storage Access Functions
120 (cffi:defcfun ("la_get_integer" ccl-la-get-integer)
121 :long (p :pointer) (i size-t)) ;; was int, not long, for first
122 (defun la-get-integer (p i)
123 (ccl-la-get-integer p i))
125 (cffi:defcfun ("la_get_double" ccl-la-get-double)
126 :double (p :pointer) (i size-t))
127 (defun la-get-double (p i)
128 (ccl-la-get-double p i))
130 (cffi:defcfun ("la_get_complex_real" ccl-la-get-complex-real)
131 :double (p :pointer) (i size-t))
132 (defun la-get-complex-real (p i)
133 (ccl-la-get-complex-real p i))
135 (cffi:defcfun ("la_get_complex_imag" ccl-la-get-complex-imag)
136 :double (p :pointer) (i size-t))
137 (defun la-get-complex-imag (p i)
138 (ccl-la-get-complex-imag p i))
140 (defun la-get-complex (p i)
141 (complex (la-get-complex-real p i) (la-get-complex-imag p i)))
143 (cffi:defcfun ("la_get_pointer" ccl-la-get-pointer)
144 :pointer (p :pointer) (i size-t))
145 (defun la-get-pointer (p i)
146 (ccl-la-get-pointer p i))
149 ;;; CFFI glue for Storage Mutation Functions
151 (cffi:defcfun ("la_put_integer" ccl-la-put-integer)
152 :void (p :pointer) (i size-t) (x :long)) ;; last was :int ?
153 (defun la-put-integer (p i x)
154 (ccl-la-put-integer p i x))
156 (cffi:defcfun ("la_put_double" ccl-la-put-double)
157 :void (p :pointer) (i size-t) (x :double))
158 (defun la-put-double (p i x)
159 (ccl-la-put-double p i (float x 1d0)))
161 (cffi:defcfun ("la_put_complex" ccl-la-put-complex)
162 :void (p :pointer) (i size-t) (x :double) (y :double))
163 (defun la-put-complex (p i x y)
164 (ccl-la-put-complex p i (float x 1d0) (float y 1d0)))
166 (cffi:defcfun ("la_put_pointer" ccl-la-put-pointer)
167 :void (p :pointer) (i size-t) (q :pointer))
168 (defun la-put-pointer (p i q)
169 (ccl-la-put-pointer p i q))
172 ;; User interface (exported)
174 (defun la-allocate (n m)
175 (let ((p (la-base-allocate n m)))
176 (if (null-ptr-p p) (error "allocation failed"))
177 (if (member p *la-allocations* :test #'ptr-eq)
178 (error "pointer is already on the list"))
179 (push p *la-allocations*)
182 (defun la-free (p)
183 (when (and (not (null-ptr-p p)) (member p *la-allocations* :test #'ptr-eq))
184 (setf *la-allocations* (delete p *la-allocations* :test #'ptr-eq))
185 (la-base-free p)))
187 (defun la-cleanup-allocations ()
188 (let ((allocs (copy-list *la-allocations*)))
189 (dolist (p allocs) (la-free p))))
191 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
192 ;;;;
193 ;;;; C Vector and Array Allocation
194 ;;;;
195 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
197 (defun la-vector(n mode) (la-allocate n (la-mode-size mode)))
198 (defun la-free-vector (v) (la-free v))
200 (defun la-matrix (n m mode)
201 (let ((matrix (la-allocate n (la-mode-size +mode-in+))))
202 (dotimes (i n)
203 (la-put-pointer matrix i (la-allocate m (la-mode-size mode))))
204 matrix))
206 (defun la-free-matrix (matrix n)
207 (dotimes (i n)
208 (la-free (la-get-pointer matrix i)))
209 (la-free matrix))
212 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
213 ;;;;
214 ;;;; C to/from Lisp Data Conversion
215 ;;;;
216 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
218 (defun la-data-to-vector (data mode)
219 (check-sequence data)
220 (let* ((n (length data))
221 (vec (la-vector n mode))
222 (d (make-next-element data)))
223 (declare (fixnum n))
224 (cond
225 ((= mode +mode-in+)
226 (dotimes (i n)
227 (declare (fixnum i))
228 (la-put-integer vec i (get-next-element d i))))
229 ((= mode +mode-re+)
230 (dotimes (i n)
231 (declare (fixnum i))
232 (la-put-double vec i (get-next-element d i))))
233 ((= mode +mode-cx+)
234 (dotimes (i n)
235 (declare (fixnum i))
236 (let ((x (get-next-element d i)))
237 (la-put-complex vec i (realpart x) (imagpart x))))))
238 vec))
240 (defun la-data-to-matrix (data mode)
241 (check-matrix data)
242 (let* ((n (num-rows data))
243 (m (num-cols data))
244 (mat (la-matrix n m mode)))
245 (declare (fixnum n m))
246 (cond
247 ((= mode +mode-in+)
248 (dotimes (i n)
249 (declare (fixnum i))
250 (let ((vec (la-get-pointer mat i)))
251 (dotimes (j m)
252 (la-put-integer vec j (aref data i j))))))
253 ((= mode +mode-re+)
254 (dotimes (i n)
255 (declare (fixnum i))
256 (let ((vec (la-get-pointer mat i)))
257 (dotimes (j m)
258 (la-put-double vec j (aref data i j))))))
259 ((= mode +mode-cx+)
260 (dotimes (i n)
261 (declare (fixnum i))
262 (let ((vec (la-get-pointer mat i)))
263 (dotimes (j m)
264 (let ((x (aref data i j)))
265 (la-put-complex vec i (realpart x) (imagpart x))))))))
266 mat))
268 (defun la-vector-to-data (vec n mode data)
269 (declare (fixnum n))
270 (check-sequence data)
271 (let ((d (make-next-element data))
272 (gf (cond
273 ((= mode +mode-in+) #'la-get-integer)
274 ((= mode +mode-re+) #'la-get-double)
275 ((= mode +mode-cx+) #'la-get-complex))))
276 (dotimes (i n)
277 (declare (fixnum i))
278 (set-next-element d i (funcall gf vec i))))
279 data)
281 (defun la-matrix-to-data (mat n m mode result)
282 (declare (fixnum n m))
283 (check-matrix result)
284 (let ((gf (cond
285 ((= mode +mode-in+) #'la-get-integer)
286 ((= mode +mode-re+) #'la-get-double)
287 ((= mode +mode-cx+) #'la-get-complex))))
288 (dotimes (i n)
289 (declare (fixnum i))
290 (let ((vec (la-get-pointer mat i)))
291 (dotimes (j m)
292 (declare (fixnum j))
293 (setf (aref result i j) (funcall gf vec j))))))
294 result)