Support for different matrix packages. Pick one and use it. BUT SBCL is only curre...
[CommonLispStat.git] / ladata.lsp
blob418c6f501a000a6238143491fec2e3b4969897c7
1 ;;; -*- mode: lisp -*-
2 ;;; Copyright (c) 2005--2007, 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 ;;;
12 ;;; Package Setup
13 ;;;
15 ;;(in-package #:lisp-stat-basics)
16 ;;(in-package :cl-user)
18 (defpackage :lisp-stat-linalg-data
19 (:use :common-lisp
20 :cffi
21 :lisp-stat-ffi-int
22 :lisp-stat-types
23 :lisp-stat-compound-data
24 :lisp-stat-matrix)
25 (:export ;; more to add
26 +mode-in+ +mode-re+ +mode-cx+ mode-of
28 la-data-mode la-allocate la-free
30 la-get-double la-put-double
31 la-put-integer
33 la-matrix la-free-matrix la-matrix-to-data la-data-to-matrix
34 la-vector la-free-vector la-vector-to-data la-data-to-vector ))
37 (in-package :lisp-stat-linalg-data)
42 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
43 ;;;
44 ;;; Data Mode Functions
45 ;;;
46 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
48 ;;;;
49 ;;;; These constants need to be redefined if IN, RE or CX in linalg.h
50 ;;;; are changed.
51 ;;;;
53 ;;; FIXME:AJR: This is how Luke got around having appropriate
54 ;;; approaches for Linear Algebra. We want to cheat and instead use
55 ;;; CLEM or MATLISP as the underlying linear algebra package.
57 (defparameter +mode-in+ 0)
58 (defparameter +mode-re+ 1)
59 (defparameter +mode-cx+ 2)
61 (defun mode-of (x)
62 (etypecase x
63 (fixnum +mode-in+)
64 (rational +mode-re+)
65 (float +mode-re+)
66 (complex +mode-cx+)))
68 (defun la-data-mode (data)
69 (let ((data (compound-data-seq data))
70 (mode 0))
71 (cond
72 ((vectorp data)
73 (let ((n (length data)))
74 (declare (fixnum n))
75 (dotimes (i n mode)
76 (declare (fixnum i))
77 (setf mode (max mode (mode-of (aref data i)))))))
78 ((consp data) (dolist (x data mode) (setf mode (max mode (mode-of x)))))
79 (t (error "bad sequence - ~s" data)))))
82 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
83 ;;;;
84 ;;;; Internal Allocation Funcitons
85 ;;;;
86 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
88 (defvar *la-allocations* nil)
90 ;;;
91 ;;; CFFI glue for... Storage Allocation Functions
92 ;;;
94 (defun null-ptr-p (p) (cffi:null-pointer-p p))
95 (defun ptr-eq (p q) (cffi:pointer-eq p q))
97 (cffi:defcfun ("la_base_allocate" ccl-la-base-allocate)
98 :pointer (n :int) (m :int))
99 (defun la-base-allocate (n m)
100 (ccl-la-base-allocate n m))
102 (cffi:defcfun ("la_base_free_alloc" ccl-la-base-free-alloc)
103 :void (p :pointer))
104 (defun la-base-free (p)
105 (ccl-la-base-free-alloc p))
107 (cffi:defcfun ("la_mode_size" ccl-la-mode-size)
108 :int (x :int))
110 (defun la-mode-size (mode)
111 (ccl-la-mode-size mode))
114 ;;; Callbacks for Internal Storage
117 (cffi:defcallback lisp-la-allocate :void ((n :long) (m :long))
118 (ccl-store-ptr (la-allocate n m)))
119 (cffi:defcfun ("register_la_allocate" register-la-allocate)
120 :void (p :pointer))
121 (register-la-allocate (cffi:callback lisp-la-allocate))
122 (cffi:defcfun ("la_allocate" la)
123 :pointer (x :int) (y :int))
125 (cffi:defcallback lisp-la-free-alloc
126 :void ((p :pointer))
127 (la-free p))
129 (cffi:defcfun ("register_la_free_alloc" register-la-free-alloc)
130 :void (p :pointer))
131 (register-la-free-alloc (cffi:callback lisp-la-free-alloc))
132 (cffi:defcfun ("la_free_alloc" lf)
133 :void (p :pointer))
137 ;;; CFFI glue for... Storage Access Functions
140 (cffi:defcfun ("la_get_integer" ccl-la-get-integer)
141 :int (p :pointer) (i :int))
142 (defun la-get-integer (p i)
143 (ccl-la-get-integer p i))
145 (cffi:defcfun ("la_get_double" ccl-la-get-double)
146 :double (p :pointer) (i :int))
147 (defun la-get-double (p i)
148 (ccl-la-get-double p i))
150 (cffi:defcfun ("la_get_complex_real" ccl-la-get-complex-real)
151 :double (p :pointer) (i :int))
152 (defun la-get-complex-real (p i)
153 (ccl-la-get-complex-real p i))
155 (cffi:defcfun ("la_get_complex_imag" ccl-la-get-complex-imag)
156 :double (p :pointer) (i :int))
157 (defun la-get-complex-imag (p i)
158 (ccl-la-get-complex-imag p i))
160 (defun la-get-complex (p i)
161 (complex (la-get-complex-real p i) (la-get-complex-imag p i)))
163 (cffi:defcfun ("la_get_pointer" ccl-la-get-pointer)
164 :pointer (p :pointer) (i :int))
165 (defun la-get-pointer (p i)
166 (ccl-la-get-pointer p i))
169 ;;; CFFI glue for Storage Mutation Functions
171 (cffi:defcfun ("la_put_integer" ccl-la-put-integer)
172 :void (p :pointer) (i :int) (x :int))
173 (defun la-put-integer (p i x)
174 (ccl-la-put-integer p i x))
176 (cffi:defcfun ("la_put_double" ccl-la-put-double)
177 :void (p :pointer) (i :int) (x :double))
178 (defun la-put-double (p i x)
179 (ccl-la-put-double p i (float x 1d0)))
181 (cffi:defcfun ("la_put_complex" ccl-la-put-complex)
182 :void (p :pointer) (i :int) (x :double) (y :double))
183 (defun la-put-complex (p i x y)
184 (ccl-la-put-complex p i (float x 1d0) (float y 1d0)))
186 (cffi:defcfun ("la_put_pointer" ccl-la-put-pointer)
187 :void (p :pointer) (i :int) (q :pointer))
188 (defun la-put-pointer (p i q)
189 (ccl-la-put-pointer p i q))
192 ;; User interface (exported)
194 (defun la-allocate (n m)
195 (let ((p (la-base-allocate n m)))
196 (if (null-ptr-p p) (error "allocation failed"))
197 (if (member p *la-allocations* :test #'ptr-eq)
198 (error "pointer is already on the list"))
199 (push p *la-allocations*)
202 (defun la-free (p)
203 (when (and (not (null-ptr-p p)) (member p *la-allocations* :test #'ptr-eq))
204 (setf *la-allocations* (delete p *la-allocations* :test #'ptr-eq))
205 (la-base-free p)))
207 (defun la-cleanup-allocations ()
208 (let ((allocs (copy-list *la-allocations*)))
209 (dolist (p allocs) (la-free p))))
211 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
212 ;;;;
213 ;;;; C Vector and Array Allocation
214 ;;;;
215 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
217 (defun la-vector(n mode) (la-allocate n (la-mode-size mode)))
218 (defun la-free-vector (v) (la-free v))
220 (defun la-matrix (n m mode)
221 (let ((matrix (la-allocate n (la-mode-size +mode-in+))))
222 (dotimes (i n)
223 (la-put-pointer matrix i (la-allocate m (la-mode-size mode))))
224 matrix))
226 (defun la-free-matrix (matrix n)
227 (dotimes (i n)
228 (la-free (la-get-pointer matrix i)))
229 (la-free matrix))
232 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
233 ;;;;
234 ;;;; C to/from Lisp Data Conversion
235 ;;;;
236 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
238 (defun la-data-to-vector (data mode)
239 (check-sequence data)
240 (let* ((n (length data))
241 (vec (la-vector n mode))
242 (d (make-next-element data)))
243 (declare (fixnum n))
244 (cond
245 ((= mode +mode-in+)
246 (dotimes (i n)
247 (declare (fixnum i))
248 (la-put-integer vec i (get-next-element d i))))
249 ((= mode +mode-re+)
250 (dotimes (i n)
251 (declare (fixnum i))
252 (la-put-double vec i (get-next-element d i))))
253 ((= mode +mode-cx+)
254 (dotimes (i n)
255 (declare (fixnum i))
256 (let ((x (get-next-element d i)))
257 (la-put-complex vec i (realpart x) (imagpart x))))))
258 vec))
260 (defun la-data-to-matrix (data mode)
261 (check-matrix data)
262 (let* ((n (num-rows data))
263 (m (num-cols data))
264 (mat (la-matrix n m mode)))
265 (declare (fixnum n m))
266 (cond
267 ((= mode +mode-in+)
268 (dotimes (i n)
269 (declare (fixnum i))
270 (let ((vec (la-get-pointer mat i)))
271 (dotimes (j m)
272 (la-put-integer vec j (aref data i j))))))
273 ((= mode +mode-re+)
274 (dotimes (i n)
275 (declare (fixnum i))
276 (let ((vec (la-get-pointer mat i)))
277 (dotimes (j m)
278 (la-put-double vec j (aref data i j))))))
279 ((= mode +mode-cx+)
280 (dotimes (i n)
281 (declare (fixnum i))
282 (let ((vec (la-get-pointer mat i)))
283 (dotimes (j m)
284 (let ((x (aref data i j)))
285 (la-put-complex vec i (realpart x) (imagpart x))))))))
286 mat))
288 (defun la-vector-to-data (vec n mode data)
289 (declare (fixnum n))
290 (check-sequence data)
291 (let ((d (make-next-element data))
292 (gf (cond
293 ((= mode +mode-in+) #'la-get-integer)
294 ((= mode +mode-re+) #'la-get-double)
295 ((= mode +mode-cx+) #'la-get-complex))))
296 (dotimes (i n)
297 (declare (fixnum i))
298 (set-next-element d i (funcall gf vec i))))
299 data)
301 (defun la-matrix-to-data (mat n m mode result)
302 (declare (fixnum n m))
303 (check-matrix result)
304 (let ((gf (cond
305 ((= mode +mode-in+) #'la-get-integer)
306 ((= mode +mode-re+) #'la-get-double)
307 ((= mode +mode-cx+) #'la-get-complex))))
308 (dotimes (i n)
309 (declare (fixnum i))
310 (let ((vec (la-get-pointer mat i)))
311 (dotimes (j m)
312 (declare (fixnum j))
313 (setf (aref result i j) (funcall gf vec j))))))
314 result)