both.lisp 13.5 KB
Newer Older
1
2
;; Functions for both vectors and matrices.
;; Liam Healy 2008-04-26 20:48:44EDT both.lisp
Liam Healy's avatar
Liam Healy committed
3
;; Time-stamp: <2010-11-25 09:31:17EST both.lisp>
4
;;
Liam M. Healy's avatar
Liam M. Healy committed
5
;; Copyright 2008, 2009, 2010 Liam M. Healy
6
7
8
9
10
11
12
13
14
15
16
17
18
19
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <http://www.gnu.org/licenses/>.
20
21
22

(in-package :gsl)

23
24
25
26
;;;;****************************************************************************
;;;; Administrative (internal use)
;;;;****************************************************************************

Liam Healy's avatar
Liam Healy committed
27
(defmfun alloc-from-block ((object vector) blockptr)
28
  ("gsl_" :category :type "_alloc_from_block")
Liam Healy's avatar
Liam Healy committed
29
  ((blockptr :pointer)
30
   (0 sizet)				; offset
31
   ((size object) sizet)		; number of elements
32
   (1 sizet))				; stride
33
34
  :definition :generic
  :c-return :pointer
35
  :export nil
36
  :documentation "Allocate memory for the GSL struct given a block pointer.")
37

Liam Healy's avatar
Liam Healy committed
38
(defmfun alloc-from-block ((object matrix) blockptr)
39
  ("gsl_" :category :type "_alloc_from_block")
Liam Healy's avatar
Liam Healy committed
40
  ((blockptr :pointer)
41
42
43
44
45
   (0 sizet)				; offset
   ((first (dimensions object)) sizet)	; number of rows
   ((second (dimensions object)) sizet)	; number of columns
   ((second (dimensions object)) sizet))	; "tda" = number of columns for now
  :definition :methods
46
47
48
49
50
51
52
  :c-return :pointer
  :export nil)

;;;;****************************************************************************
;;;; Bulk operations
;;;;****************************************************************************

53
54
55
56
(defmfun set-all ((object both) value)
  ("gsl_" :category :type "_set_all")
  (((mpointer object) :pointer) (value :element-c-type))
  :definition :generic
Liam Healy's avatar
Liam Healy committed
57
  :element-types #+fsbv t #-fsbv :no-complex
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
  :outputs (object)
  :c-return :void
  :documentation "Set all elements to the value.")

(defmfun set-zero ((object both))
  ("gsl_"  :category :type "_set_zero")
  (((mpointer object) :pointer))
  :definition :generic
  :inputs (object)
  :outputs (object)
  :c-return :void
  :documentation "Set all elements to 0.")

(defmfun swap ((a both) (b both))
  ("gsl_" :category :type "_swap")
  (((mpointer a) :pointer) ((mpointer b) :pointer))
  :definition :generic
  :inputs (a b)
  :outputs (a b)
  :documentation			; FDL
  "Exchange the elements of a and b
   by copying.  The two must have the same dimensions.")

81
82
83
84
85
;;;;****************************************************************************
;;;; Array elements; used in callbacks scalarsp=T only
;;;;****************************************************************************
;;; Normal foreign array access is with grid:gref, but in order to
;;; avoid the overhead of instantiating a foreign-array object to
86
87
;;; access components, we use these macros which expand to gsl_*_get
;;; and gsl_*_set.
88

Liam Healy's avatar
Liam Healy committed
89
90
91
92
93
94
95
96
;;; It's unlikely users will need maref directly, because when using
;;; the GSL library, we use the grid functions to get elements.  The
;;; main use in the GSL library is user-defined callbacks (e.g. in
;;; solve-minimize-fit) when scalarsp = T is specified; in that case
;;; the generated wrapper function includes maref.  However, if
;;; another foreign library which used the GSL library provided an
;;; mpointer, the maref macro saves the time of creating an object
;;; around it.
97
98
99
(export 'maref)

(defun access-value-int (mpointer class-name value indices)
Liam Healy's avatar
Liam Healy committed
100
  "Create a form to access the GSL array value from the mpointer.  If value is not nil,
Liam M. Healy's avatar
Liam M. Healy committed
101
   set the value; otherwise, get the value."
102
  ;; (access-value-int 'ptr 'grid:vector-unsigned-byte-16 45 '(3))
Liam M. Healy's avatar
Liam M. Healy committed
103
  ;; (FOREIGN-FUNCALL "gsl_vector_ushort_set" :POINTER PTR SIZET 3 :UNSIGNED-SHORT 45 :VOID)
104
  ;; (access-value-int 'ptr 'grid:matrix-unsigned-byte-16 nil '(45 3))
Liam M. Healy's avatar
Liam M. Healy committed
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
  ;; (FOREIGN-FUNCALL "gsl_matrix_ushort_get" :POINTER PTR SIZET 45 SIZET 3 :UNSIGNED-SHORT)
  (let ((element-type (grid::farray-element-type class-name))
	(matrixp (subtypep class-name 'grid:matrix)))
    (unless (or (and matrixp (eql 2 (length indices)))
		(and (not matrixp) (eql 1 (length indices))))
      (error "The number of indices, ~a, is not correct for an array of class ~a"
	     (length indices)
	     class-name))
    `(cffi:foreign-funcall
      ,(actual-gsl-function-name
	`("gsl_" :category :type ,(if value "_set" "_get"))
	(if matrixp 'matrix 'vector)
	element-type)
      :pointer ,mpointer
      sizet ,(first indices)
      ,@(when matrixp (list 'sizet (second indices)))
      ,(grid:cl-cffi element-type)
      ,@(when value `(,value :void)))))

124
(defmacro maref (mpointer class-name &rest indices)
Liam Healy's avatar
Liam Healy committed
125
126
127
  "Get or set (setf maref) the array element from the GSL mpointer.
   The class-name is the specific subclass name of
   grid:foreign-array."
128
  (access-value-int mpointer class-name nil indices))
129

130
131
132
133
134
135
136
137
138
139
140
141
142
143
(defmacro set-maref (mpointer class-name &rest indices-value)
  (alexandria:once-only ((value (alexandria:lastcar indices-value)))
    `(progn
       ,(access-value-int
	 mpointer
	 class-name
	 value
	 (butlast indices-value))
       ,value)))

(defsetf maref set-maref)

;;; (maref ptr grid:vector-double-float 3)
;;; (setf (maref ptr grid:vector-double-float 3) 45.0d0)
144

145
;;;;****************************************************************************
146
;;;; Elementwise arithmetic operations overwriting an array
147
148
;;;;****************************************************************************

149
150
151
152
153
;;; These are overwriting elementwise arithmetic operations.  Where
;;; there are two array arguments, the operation acts on both arrays
;;; on corresponding elements.  The result is placed in the first
;;; array, overwriting the original contents.

Liam Healy's avatar
Liam Healy committed
154
155
156
;;; Errors in GSL:
;;; 1) complex operations in older versions of GSL
;;; https://savannah.gnu.org/bugs/index.php?22478
157
;;; Fixed in 1.12, can change the :no-complex spec.
Liam Healy's avatar
Liam Healy committed
158

159
(defmfun elt+ ((a both) (b both))
160
  ("gsl_" :category :type "_add")
161
  (((mpointer a) :pointer) ((mpointer b) :pointer))
162
  :definition :generic
163
  :element-types #.(if (have-at-least-gsl-version '(1 12)) t :no-complex)
164
165
166
167
168
169
170
  :inputs (a b)
  :outputs (a)
  :return (a)
  :documentation			; FDL
  "Add the elements of b to the elements of vector a
   The two must have the same dimensions.")

171
(defmfun elt+ ((a both) (x float))
172
173
174
  ("gsl_" :category :type "_add_constant")
  (((mpointer a) :pointer) (x :double))
  :definition :methods
175
  :element-types :no-complex
176
177
178
179
180
181
  :inputs (a)
  :outputs (a)
  :return (a)
  :documentation			; FDL
  "Add the scalar double-float x to all the elements of array a.")

182
183
184
185
186
187
188
189
190
191
192
193
194
#+fsbv
(defmfun elt+ ((a both) (x complex))
  ("gsl_" :category :type "_add_constant")
  (((mpointer a) :pointer) (x :element-c-type))
  :definition :methods
  :element-types :complex
  :inputs (a)
  :outputs (a)
  :return (a)
  :gsl-version (1 12)
  :documentation			; FDL
  "Add the scalar complex x to all the elements of array a.")

195
(defmethod elt+ ((x float) (a foreign-array))
196
  (elt+ a x))
197
  
198
(defmfun elt- ((a both) (b both))
199
200
201
  ("gsl_" :category :type "_sub")
  (((mpointer a) :pointer) ((mpointer b) :pointer))
  :definition :generic
202
  :element-types #.(if (have-at-least-gsl-version '(1 12)) t :no-complex)
203
204
205
206
207
208
  :inputs (a b)
  :outputs (a)
  :return (a)
  :documentation			; FDL
  "Subtract the elements of b from the elements of a.
   The two must have the same dimensions.")
209

210
(defmethod elt- ((a foreign-array) (x float))
211
212
  (elt+ a (- x)))

213
(defmfun elt* ((a vector) (b vector))
214
215
216
  ("gsl_" :category :type "_mul")
  (((mpointer a) :pointer) ((mpointer b) :pointer))
  :definition :generic
217
  :element-types #.(if (have-at-least-gsl-version '(1 12)) t :no-complex)
218
219
220
221
222
223
224
  :inputs (a b)
  :outputs (a)
  :return (a)
  :documentation			; FDL
  "Multiply the elements of a by the elements of b.
   The two must have the same dimensions.")

225
(defmfun elt* ((a matrix) (b matrix))
226
227
  ("gsl_" :category :type "_mul_elements")
  (((mpointer a) :pointer) ((mpointer b) :pointer))
228
  :definition :methods
Liam Healy's avatar
Liam Healy committed
229
  :element-types :no-complex
230
231
232
233
  :inputs (a b)
  :outputs (a)
  :return (a))

234
(defmfun elt/ ((a vector) (b vector))
235
236
237
  ("gsl_" :category :type "_div")
  (((mpointer a) :pointer) ((mpointer b) :pointer))
  :definition :generic
238
  :element-types #.(if (have-at-least-gsl-version '(1 12)) t :no-complex)
239
240
241
242
243
244
245
  :inputs (a b)
  :outputs (a)
  :return (a)
  :documentation			; FDL
  "Divide the elements of a by the elements of b.
   The two must have the same dimensions.")

246
(defmfun elt/ ((a matrix) (b matrix))
247
248
  ("gsl_" :category :type "_div_elements")
  (((mpointer a) :pointer) ((mpointer b) :pointer))
249
  :definition :methods
Liam Healy's avatar
Liam Healy committed
250
  :element-types :no-complex
251
252
253
254
  :inputs (a b)
  :outputs (a)
  :return (a))

255
(defmethod elt/ ((a foreign-array) (x number))
256
257
  (elt* a (/ x)))

258
(defmfun elt* ((a both) (x float))
259
  ("gsl_" :category :type "_scale")
260
  (((mpointer a) :pointer) (x :double))
261
  :definition :methods
262
  :element-types :no-complex
263
264
265
266
  :inputs (a)
  :outputs (a)
  :return (a)
  :documentation			; FDL
267
  "Multiply the elements of a by the scalar double-float factor x.")
268

269
270
271
272
273
274
275
276
277
278
279
280
281
#+fsbv
(defmfun elt* ((a both) (x complex))
  ("gsl_" :category :type "_scale")
  (((mpointer a) :pointer) (x :element-c-type))
  :definition :methods
  :element-types :complex
  :inputs (a)
  :outputs (a)
  :return (a)
  :gsl-version (1 12)
  :documentation			; FDL
  "Multiply the elements of a by the scalar complex factor x.")

282
(defmethod elt* ((x float) (a foreign-array))
283
  (elt* a x))
284
285
286
287
288
289
290
291
292

;;;;****************************************************************************
;;;; Maximum and minimum elements
;;;;****************************************************************************

(defmfun mmax ((a both))
  ("gsl_" :category :type "_max")
  (((mpointer a) :pointer))
  :definition :generic
293
  :element-types :no-complex
294
295
296
297
298
299
300
301
302
  :inputs (a)
  :c-return :element-c-type
  :documentation			; FDL
  "The maximum value in a.")

(defmfun mmin ((a both))
  ("gsl_" :category :type "_min")
  (((mpointer a) :pointer))
  :definition :generic
303
  :element-types :no-complex
304
305
306
307
308
309
310
  :inputs (a)
  :c-return :element-c-type
  :documentation			; FDL
  "The minimum value in a.")

(defmfun minmax ((a both))
  ("gsl_" :category :type "_minmax")
311
312
313
  (((mpointer a) :pointer)
   (min (:pointer :element-c-type))
   (max (:pointer :element-c-type)))
314
  :definition :generic
315
  :element-types :no-complex
316
317
318
319
320
321
322
323
324
  :inputs (a)
  :c-return :void
  :documentation			; FDL
  "The minimum and maximum values in a.")

(defmfun min-index ((a vector))
  ("gsl_" :category :type "_min_index")
  (((mpointer a) :pointer))
  :definition :generic
325
  :element-types :no-complex
326
  :inputs (a)
327
  :c-return sizet
328
329
330
331
332
333
  :documentation			; FDL
  "The index of the minimum value in a.  When there are several
  equal minimum elements, then the lowest index is returned.")

(defmfun min-index ((a matrix))
  ("gsl_" :category :type "_min_index")
334
  (((mpointer a) :pointer) (imin (:pointer sizet)) (jmin (:pointer sizet)))
335
  :definition :methods
336
  :element-types :no-complex
337
338
339
340
341
342
343
  :inputs (a)
  :c-return :void)

(defmfun max-index ((a vector))
  ("gsl_" :category :type "_max_index")
  (((mpointer a) :pointer))
  :definition :generic
344
  :element-types :no-complex
345
  :inputs (a)
346
  :c-return sizet
347
348
349
350
351
352
  :documentation			; FDL
  "The index of the maximum value in a.  When there are several
  equal maximum elements, then the lowest index is returned.")

(defmfun max-index ((a matrix))
  ("gsl_" :category :type "_max_index")
353
  (((mpointer a) :pointer) (imin (:pointer sizet)) (jmin (:pointer sizet)))
354
  :definition :methods
355
  :element-types :no-complex
356
357
358
359
360
  :inputs (a)
  :c-return :void)

(defmfun minmax-index ((a vector))
  ("gsl_" :category :type "_minmax_index")
361
  (((mpointer a) :pointer) (imin (:pointer sizet)) (jmin (:pointer sizet)))
362
  :definition :generic
363
  :element-types :no-complex
364
365
366
367
368
369
370
371
372
373
  :inputs (a)
  :c-return :void
  :documentation			; FDL
  "The indices of the minimum and maximum values in a.
  When there are several equal minimum elements then the lowest index is
  returned.  Returned indices are minimum, maximum; for matrices
  imin, jmin, imax, jmax.")

(defmfun minmax-index ((a matrix))
  ("gsl_" :category :type "_minmax_index")
374
375
376
  (((mpointer a) :pointer)
   (imin (:pointer sizet)) (jmin (:pointer sizet))
   (imax (:pointer sizet)) (jmax (:pointer sizet)))
377
  :definition :methods
378
  :element-types :no-complex
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
  :inputs (a)
  :c-return :void)

;;;;****************************************************************************
;;;; Properties
;;;;****************************************************************************

(defmfun mzerop ((a both))
  ("gsl_" :category :type "_isnull")
  (((mpointer a) :pointer))
  :definition :generic
  :inputs (a)
  :c-return :boolean
  :documentation			; FDL
  "All elements of a are zero.")
Liam Healy's avatar
Liam Healy committed
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423

(defmfun mplusp ((a both))
  ("gsl_" :category :type "_ispos")
  (((mpointer a) :pointer))
  :definition :generic
  :gsl-version (1 9)
  :inputs (a)
  :c-return :boolean
  :documentation			; FDL
  "All elements of a are positive.")

(defmfun mminusp ((a both))
  ("gsl_" :category :type "_isneg")
  (((mpointer a) :pointer))
  :definition :generic
  :gsl-version (1 9)
  :inputs (a)
  :c-return :boolean
  :documentation			; FDL
  "All elements of a are negative.")

(defmfun non-negative-p ((a both))
  ("gsl_" :category :type "_isnonneg")
  (((mpointer a) :pointer))
  :definition :generic
  :gsl-version (1 10)
  :inputs (a)
  :c-return :boolean
  :documentation			; FDL
  "All elements of a are non-negative.")