src/numerics/ladata.lsp

   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.
   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 ;;;
  12 ;;; Package Setup
  13 ;;;
  14
  15 ;;(in-package #:lisp-stat-basics)
  16 ;;(in-package :cl-user)
  17
  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
  27
  28       la-data-mode la-allocate la-free
  29
  30       la-get-double la-put-double
  31       la-put-integer
  32
  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 ))
  35
  36
  37 (in-package :lisp-stat-linalg-data)
  38
  39 #+openmcl
  40 (defctype size-t :unsigned-long)
  41 #+sbcl
  42 (defctype size-t :unsigned-int)
  43
  44 ;; Should we do the same with int's and long's?   There is some
  45 ;; evidence that it might be useful?
  46
  47
  48 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  49 ;;;
  50 ;;;                        Data Mode Functions
  51 ;;;
  52 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  53
  54 ;;;;
  55 ;;;; These constants need to be redefined if IN, RE or CX in linalg.h
  56 ;;;; are changed.
  57 ;;;;
  58
  59 ;;; FIXME:AJR: This is how Luke got around having appropriate
  60 ;;; approaches for Linear Algebra.  We want to cheat and instead use
  61 ;;; CLEM or MATLISP as the underlying linear algebra package.
  62
  63 (defparameter +mode-in+ 0)
  64 (defparameter +mode-re+ 1)
  65 (defparameter +mode-cx+ 2)
  66
  67 (defun mode-of (x)
  68   (etypecase x
  69              (fixnum +mode-in+)
  70              (rational +mode-re+)
  71              (float +mode-re+)
  72              (complex +mode-cx+)))
  73
  74 (defun la-data-mode (data)
  75   (let ((data (compound-data-seq data))
  76         (mode 0))
  77     (cond
  78      ((vectorp data)
  79       (let ((n (length data)))
  80         (declare (fixnum n))
  81         (dotimes (i n mode)
  82                  (declare (fixnum i))
  83                  (setf mode (max mode (mode-of (aref data i)))))))
  84      ((consp data) (dolist (x data mode) (setf mode (max mode (mode-of x)))))
  85      (t (error "bad sequence - ~s" data)))))
  86
  87
  88 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  89 ;;;;
  90 ;;;;                        Internal Allocation Funcitons
  91 ;;;;
  92 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  93
  94 (defvar *la-allocations* nil)
  95
  96 ;;;
  97 ;;; CFFI glue for...     Storage Allocation Functions
  98 ;;;
  99
 100 (defun null-ptr-p (p) (cffi:null-pointer-p p))
 101 (defun ptr-eq (p q) (cffi:pointer-eq p q))
 102
 103 (cffi:defcfun ("la_base_allocate" ccl-la-base-allocate)
 104     :pointer (n size-t) (m size-t))
 105 (defun la-base-allocate (n m)
 106   (ccl-la-base-allocate n m))
 107
 108 (cffi:defcfun ("la_base_free_alloc" ccl-la-base-free-alloc)
 109     :void (p :pointer))
 110 (defun la-base-free (p)
 111   (ccl-la-base-free-alloc p))
 112
 113 (cffi:defcfun ("la_mode_size" ccl-la-mode-size)
 114     size-t (x size-t))
 115
 116 (defun la-mode-size (mode)
 117   (ccl-la-mode-size mode))
 118
 119 ;;;
 120 ;;;             Callbacks for Internal Storage
 121 ;;;
 122
 123 (cffi:defcallback lisp-la-allocate :void ((n size-t) (m size-t))
 124                   (ccl-store-ptr (la-allocate n m)))
 125 (cffi:defcfun ("register_la_allocate" register-la-allocate)
 126     :void (p :pointer))
 127 (register-la-allocate (cffi:callback lisp-la-allocate))
 128 (cffi:defcfun ("la_allocate" la)
 129     :pointer (x size-t) (y size-t))
 130
 131 (cffi:defcallback lisp-la-free-alloc
 132     :void ((p :pointer))
 133   (la-free p))
 134
 135 (cffi:defcfun ("register_la_free_alloc" register-la-free-alloc)
 136     :void (p :pointer))
 137 (register-la-free-alloc (cffi:callback lisp-la-free-alloc))
 138 (cffi:defcfun ("la_free_alloc" lf)
 139     :void (p :pointer))
 140
 141
 142
 143 ;;; CFFI glue for...         Storage Access Functions
 144
 145
 146 (cffi:defcfun ("la_get_integer" ccl-la-get-integer)
 147     :long (p :pointer) (i size-t))   ;; was int, not long, for first
 148 (defun la-get-integer (p i)
 149   (ccl-la-get-integer p i))
 150
 151 (cffi:defcfun ("la_get_double" ccl-la-get-double)
 152     :double (p :pointer) (i size-t))
 153 (defun la-get-double (p i)
 154   (ccl-la-get-double p i))
 155
 156 (cffi:defcfun ("la_get_complex_real" ccl-la-get-complex-real)
 157     :double (p :pointer) (i size-t))
 158 (defun la-get-complex-real (p i)
 159   (ccl-la-get-complex-real p i))
 160
 161 (cffi:defcfun ("la_get_complex_imag" ccl-la-get-complex-imag)
 162     :double (p :pointer) (i size-t))
 163 (defun la-get-complex-imag (p i)
 164   (ccl-la-get-complex-imag p i))
 165
 166 (defun la-get-complex (p i)
 167   (complex (la-get-complex-real p i) (la-get-complex-imag p i)))
 168
 169 (cffi:defcfun ("la_get_pointer" ccl-la-get-pointer)
 170     :pointer (p :pointer) (i size-t))
 171 (defun la-get-pointer (p i)
 172   (ccl-la-get-pointer p i))
 173
 174
 175 ;;; CFFI glue for     Storage Mutation Functions
 176
 177 (cffi:defcfun ("la_put_integer" ccl-la-put-integer)
 178     :void (p :pointer) (i size-t) (x :long))  ;; last was :int ?
 179 (defun la-put-integer (p i x)
 180   (ccl-la-put-integer p i x))
 181
 182 (cffi:defcfun ("la_put_double" ccl-la-put-double)
 183     :void (p :pointer) (i size-t) (x :double))
 184 (defun la-put-double (p i x)
 185   (ccl-la-put-double p i (float x 1d0)))
 186
 187 (cffi:defcfun ("la_put_complex" ccl-la-put-complex)
 188     :void (p :pointer) (i size-t) (x :double) (y :double))
 189 (defun la-put-complex (p i x y)
 190   (ccl-la-put-complex p i (float x 1d0) (float y 1d0)))
 191
 192 (cffi:defcfun ("la_put_pointer" ccl-la-put-pointer)
 193     :void (p :pointer) (i size-t) (q :pointer))
 194 (defun la-put-pointer (p i q)
 195   (ccl-la-put-pointer p i q))
 196
 197
 198 ;; User interface (exported)
 199
 200 (defun la-allocate (n m)
 201   (let ((p (la-base-allocate n m)))
 202     (if (null-ptr-p p) (error "allocation failed"))
 203     (if (member p *la-allocations* :test #'ptr-eq)
 204       (error "pointer is already on the list"))
 205     (push p *la-allocations*)
 206     p))
 207
 208 (defun la-free (p)
 209   (when (and (not (null-ptr-p p)) (member p *la-allocations* :test #'ptr-eq))
 210         (setf *la-allocations* (delete p *la-allocations* :test #'ptr-eq))
 211         (la-base-free p)))
 212
 213 (defun la-cleanup-allocations ()
 214   (let ((allocs (copy-list *la-allocations*)))
 215     (dolist (p allocs) (la-free p))))
 216
 217 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 218 ;;;;
 219 ;;;;                       C Vector and Array Allocation
 220 ;;;;
 221 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 222
 223 (defun la-vector(n mode) (la-allocate n (la-mode-size mode)))
 224 (defun la-free-vector (v) (la-free v))
 225
 226 (defun la-matrix (n m mode)
 227   (let ((matrix (la-allocate n (la-mode-size +mode-in+))))
 228     (dotimes (i n)
 229       (la-put-pointer matrix i (la-allocate m (la-mode-size mode))))
 230     matrix))
 231
 232 (defun la-free-matrix (matrix n)
 233   (dotimes (i n)
 234     (la-free (la-get-pointer matrix i)))
 235   (la-free matrix))
 236
 237
 238 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 239 ;;;;
 240 ;;;;                   C to/from Lisp Data Conversion
 241 ;;;;
 242 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 243
 244 (defun la-data-to-vector (data mode)
 245   (check-sequence data)
 246   (let* ((n (length data))
 247          (vec (la-vector n mode))
 248          (d (make-next-element data)))
 249     (declare (fixnum n))
 250     (cond
 251      ((= mode +mode-in+)
 252       (dotimes (i n)
 253         (declare (fixnum i))
 254         (la-put-integer vec i (get-next-element d i))))
 255      ((= mode +mode-re+)
 256       (dotimes (i n)
 257         (declare (fixnum i))
 258         (la-put-double vec i (get-next-element d i))))
 259       ((= mode +mode-cx+)
 260        (dotimes (i n)
 261          (declare (fixnum i))
 262          (let ((x (get-next-element d i)))
 263            (la-put-complex vec i (realpart x) (imagpart x))))))
 264     vec))
 265
 266 (defun la-data-to-matrix (data mode)
 267   (check-matrix data)
 268   (let* ((n (num-rows data))
 269          (m (num-cols data))
 270          (mat (la-matrix n m mode)))
 271     (declare (fixnum n m))
 272     (cond
 273      ((= mode +mode-in+)
 274       (dotimes (i n)
 275         (declare (fixnum i))
 276         (let ((vec (la-get-pointer mat i)))
 277           (dotimes (j m)
 278             (la-put-integer vec j (aref data i j))))))
 279       ((= mode +mode-re+)
 280        (dotimes (i n)
 281          (declare (fixnum i))
 282          (let ((vec (la-get-pointer mat i)))
 283            (dotimes (j m)
 284              (la-put-double vec j (aref data i j))))))
 285       ((= mode +mode-cx+)
 286        (dotimes (i n)
 287          (declare (fixnum i))
 288          (let ((vec (la-get-pointer mat i)))
 289            (dotimes (j m)
 290              (let ((x (aref data i j)))
 291                (la-put-complex vec i (realpart x) (imagpart x))))))))
 292     mat))
 293
 294 (defun la-vector-to-data (vec n mode data)
 295   (declare (fixnum n))
 296   (check-sequence data)
 297   (let ((d (make-next-element data))
 298         (gf (cond
 299               ((= mode +mode-in+) #'la-get-integer)
 300               ((= mode +mode-re+) #'la-get-double)
 301               ((= mode +mode-cx+) #'la-get-complex))))
 302     (dotimes (i n)
 303       (declare (fixnum i))
 304       (set-next-element d i (funcall gf vec i))))
 305   data)
 306
 307 (defun la-matrix-to-data (mat n m mode result)
 308   (declare (fixnum n m))
 309   (check-matrix result)
 310   (let ((gf (cond
 311               ((= mode +mode-in+) #'la-get-integer)
 312               ((= mode +mode-re+) #'la-get-double)
 313               ((= mode +mode-cx+) #'la-get-complex))))
 314     (dotimes (i n)
 315       (declare (fixnum i))
 316       (let ((vec (la-get-pointer mat i)))
 317         (dotimes (j m)
 318           (declare (fixnum j))
 319           (setf (aref result i j) (funcall gf vec j))))))
 320   result)
 321