src/numerics/ladata.lsp

   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.
   5
   6 ;;;; ladata -- Data handling functions for linear algebra interface
   7 ;;;;
   8 ;;;; Copyright (c) 1991, by Luke Tierney. Permission is granted for
   9 ;;;; unrestricted use.
  10
  11 (in-package :lisp-stat-linalg-data)
  12
  13 #+openmcl
  14 (defctype size-t :unsigned-long)
  15 #+sbcl
  16 (defctype size-t :unsigned-int)
  17
  18 ;; Should we do the same with int's and long's?   There is some
  19 ;; evidence that it might be useful?
  20
  21
  22 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  23 ;;;
  24 ;;;                        Data Mode Functions
  25 ;;;
  26 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  27
  28 ;;;;
  29 ;;;; These constants need to be redefined if IN, RE or CX in linalg.h
  30 ;;;; are changed.
  31 ;;;;
  32
  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.
  36
  37 (defparameter +mode-in+ 0)
  38 (defparameter +mode-re+ 1)
  39 (defparameter +mode-cx+ 2)
  40
  41 (defun mode-of (x)
  42   (etypecase x
  43     (fixnum +mode-in+)
  44     (rational +mode-re+)
  45     (float +mode-re+)
  46     (complex +mode-cx+)))
  47
  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)))))
  60
  61
  62 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  63 ;;;;
  64 ;;;;                        Internal Allocation Funcitons
  65 ;;;;
  66 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  67
  68 (defvar *la-allocations* nil)
  69
  70 ;;;
  71 ;;; CFFI glue for...     Storage Allocation Functions
  72 ;;;
  73
  74 (defun null-ptr-p (p) (cffi:null-pointer-p p))
  75 (defun ptr-eq (p q) (cffi:pointer-eq p q))
  76
  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))
  81
  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))
  86
  87 (cffi:defcfun ("la_mode_size" ccl-la-mode-size)
  88     size-t (x size-t))
  89
  90 (defun la-mode-size (mode)
  91   (ccl-la-mode-size mode))
  92
  93 ;;;
  94 ;;;             Callbacks for Internal Storage
  95 ;;;
  96
  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))
 104
 105 (cffi:defcallback lisp-la-free-alloc
 106     :void ((p :pointer))
 107   (la-free p))
 108
 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))
 114
 115
 116
 117 ;;; CFFI glue for...         Storage Access Functions
 118
 119
 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))
 124
 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))
 129
 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))
 134
 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))
 139
 140 (defun la-get-complex (p i)
 141   (complex (la-get-complex-real p i) (la-get-complex-imag p i)))
 142
 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))
 147
 148
 149 ;;; CFFI glue for     Storage Mutation Functions
 150
 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))
 155
 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)))
 160
 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)))
 165
 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))
 170
 171
 172 ;; User interface (exported)
 173
 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*)
 180     p))
 181
 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)))
 186
 187 (defun la-cleanup-allocations ()
 188   (let ((allocs (copy-list *la-allocations*)))
 189     (dolist (p allocs) (la-free p))))
 190
 191 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 192 ;;;;
 193 ;;;;                       C Vector and Array Allocation
 194 ;;;;
 195 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 196
 197 (defun la-vector(n mode) (la-allocate n (la-mode-size mode)))
 198 (defun la-free-vector (v) (la-free v))
 199
 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))
 205
 206 (defun la-free-matrix (matrix n)
 207   (dotimes (i n)
 208     (la-free (la-get-pointer matrix i)))
 209   (la-free matrix))
 210
 211
 212 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 213 ;;;;
 214 ;;;;                   C to/from Lisp Data Conversion
 215 ;;;;
 216 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 217
 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))
 239
 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))
 267
 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)
 280
 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)
 295