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