Skip to content
Snippets Groups Projects
Commit 386b4440 authored by lhealy's avatar lhealy Committed by Liam Healy
Browse files

Remove obsolete files for data.

git-svn-id: svn+ssh://common-lisp.net/project/gsll/svn/branches/ffa@62 af03a46b-e846-0410-96e5-d5653c316fd0
parent 5eaf5234
No related branches found
No related tags found
No related merge requests found
;; Blocks of data
;; Liam Healy, Mon Mar 27 2006 - 12:28
;; Time-stamp: <2008-03-09 18:54:10EDT block.lisp>
;; $Id$
(in-package :gsl)
;;; Block definition
(cffi:defcstruct block
(size size)
(data :pointer))
(add-data-class block double-float block-double-float block "block")
(add-data-class block single-float block-single-float block "block")
(add-data-class block fixnum block-fixnum block "block")
(add-data-class block complex block-complex block "block")
(defdata block double-float)
(defdata block single-float)
(defdata block fixnum)
(defdata block complex)
;; Expansions for GSL data objects
;; Liam Healy 2008-03-23 14:22:01EDT data-expansions.lisp
;; Time-stamp: <2008-03-23 18:06:22EDT data-expansions.lisp>
;; $Id$
(in-package :gsl)
;;; Replace the data-go macro from data.lisp
(defmacro data-go (type matrixp)
"Define the letm function for data types."
)
;; Using GSL bulk data (vectors, matrices, etc.) storage.
;; Liam Healy, Sun Mar 26 2006 - 16:32
;; Time-stamp: <2008-03-30 14:16:25EDT data.lisp>
;; $Id$
(in-package :gsl)
;;; To do:
;;; - create a CL object of class gsl-data with the right GSL
;;; pointer just from a raw CL object
;;; - recreate GSL object should be possible to recreate the C object
;;; (need this?)
;;; - master list of objects, for manual memory management?
;;;;****************************************************************************
;;;; Class gsl-data and generic functions
;;;;****************************************************************************
(defclass gsl-data ()
((pointer :initarg :pointer :accessor pointer
:documentation "A C pointer to the GSL representation of the data.")
(storage-size :initarg :storage-size :reader storage-size)
(data :accessor data-cache
:documentation "The Lisp object corresponding to the GSL data.")
(cl-invalid
:initform t :accessor cl-invalid
:documentation
"An indication of whether the Lisp object (slot 'data) agrees with the
GSL C data. If NIL, they agree. If T, they disagree in an unspecified
way. If a list of index sets, those indices disagree and the remainder
are correct."))
(:documentation
"A superclass for all GSL data storage structures, such as vector, matrix,
etc."))
(defmethod print-object ((object gsl-data) stream)
(print-data-object object *print-array* stream))
(defparameter *print-data-contents* t)
(defun print-data-object (object contents stream)
"Print the data object to the stream, possibly showing contents."
(print-unreadable-object (object stream :type t :identity t)
(when (and contents *print-data-contents*)
(princ (data object) stream))))
(defgeneric gsl-array (object)
(:documentation "A pointer to the GSL array with the data contents."))
(defun dim0 (object)
"The first dimension of the object."
(first (storage-size object)))
(defun dim1 (object)
"The second dimension of the object."
(second (storage-size object)))
(defgeneric cl-elt-type (object)
(:documentation "The CL type of an element."))
;;; Accessing elements
(export 'maref)
(defgeneric maref (object &rest indices)
(:documentation "An element of the data."))
(defgeneric (setf maref) (value object &rest indices)
(:method :after (value (object gsl-data) &rest indices)
(push indices (cl-invalid object)))
(:documentation "Set an element of the data."))
(export '(write-binary read-binary write-formatted read-formatted))
(defgeneric write-binary (object stream)
(:documentation "Write the binary GSL data."))
(defgeneric read-binary (object stream)
(:documentation "Read the binary GSL data."))
(defgeneric write-formatted (object stream format)
(:documentation "Write the formatted GSL data."))
(defgeneric read-formatted (object stream format)
(:documentation "Read the formatted GSL data."))
;;;;****************************************************************************
;;;; Macro defdata to define class etc.
;;;;****************************************************************************
(defun assign-pointer (object pointer)
"Check that a GSL data pointer is not null, then assign it to the object."
(check-null-pointer
pointer
:ENOMEM
(format nil "for ~a."
(with-output-to-string (stream)
(print-data-object object nil stream))))
(setf (pointer object) pointer))
;;; (args (loop for i below dimensions collect (intern (format nil "I~d" i))))
;;; (mapcar (lambda (v) `(,v size)) args)
(defmacro data-go (type matrixp)
"Define the letm function for data types."
(if matrixp ; Matrix (two indices)
`(defgo ,type (size-or-initial &optional size-or-zero zero)
(cond ((and (numberp size-or-initial) (numberp size-or-zero))
;; Array dimensions are given as literal numbers for matrix
(list
`(make-data ',',type ,zero ,size-or-initial ,size-or-zero)
'free))
((typep size-or-initial 'array)
(let ((argsymb (gensym "ARG")))
(list
`(apply #'make-data
',',type ,size-or-zero (array-dimensions ,argsymb))
'free
(lambda (symb)
`(setf (data ,symb) ,argsymb))
(lambda () (list argsymb size-or-initial)))))
(t
;; Determine at runtime whether first arg is initial or dimension
(let ((argsymb (gensym "ARG")))
(list
`(apply #'make-data ',',type nil
(if (and (numberp ,argsymb) (numberp ,size-or-zero))
(list ,argsymb ,size-or-zero)
(array-dimensions ,argsymb)))
'free
(lambda (symb)
`(unless (and (numberp ,argsymb) (numberp ,size-or-zero))
(setf (data ,symb) ,argsymb)))
(lambda () (list argsymb size-or-initial)))))))
;; Vector or other one-index object
`(defgo ,type (size-or-initial &optional zero)
(typecase size-or-initial
(number
;; Size is given as literal number
(list
`(make-data ',',type ,zero ,size-or-initial)
'free))
(vector
;; Initial value supplied as literal CL vector)
(let ((argsymb (gensym "ARG")))
(list
`(make-data ',',type ,zero
(length ,argsymb))
'free
(lambda (symb) `(setf (data ,symb) ,argsymb))
(lambda () (list argsymb size-or-initial)))))
(t
;; Determine at runtime whether first arg is initial or size
(let ((argsymb (gensym "ARG")))
(list
`(make-data ',',type ,zero
(if (numberp ,argsymb) ,argsymb (length ,argsymb)))
'free
(lambda (symb) `(unless (numberp ,argsymb) (setf (data ,symb) ,argsymb)))
(lambda () (list argsymb size-or-initial)))))))))
(defparameter *data-class-name* nil
"A list classes, each consisting of a list
superclass, CL element type, class, GSL splice name.")
(defmacro add-data-class (category element-type class superclass GSL-string)
`(eval-when (:compile-toplevel :load-toplevel :execute)
(pushnew (list ',category ',element-type ',class ',superclass ',GSL-string)
*data-class-name*
:test #'equal)))
(defun data-type-lookup (category element-type)
(find (list category element-type)
*data-class-name*
:key (lambda (l) (subseq l 0 2))
:test #'equal))
(defun data-class-name (category element-type)
(third (data-type-lookup category element-type)))
(defun data-superclass-name (category element-type)
(fourth (data-type-lookup category element-type)))
(defun data-gsl-string (category element-type)
(fifth (data-type-lookup category element-type)))
(defmacro defdata (category cl-elt-type &optional (dimensions 1) splice-name)
"For the type named in the string,
define the allocator (gsl-*-alloc), zero allocator (gsl-*-calloc),
freeing (gsl-*-free), binary writing (binary-*-write) and
reading (binary-*-read), formatted writing (write-*-formatted)
and reading (read-*-formatted) functions."
(flet ((gsl-name (function-name)
(format nil "gsl_~a~a_~a"
(data-gsl-string category cl-elt-type)
(or splice-name
(lookup-splice-name cl-elt-type)) function-name)))
(let* ((cargs (loop for i below dimensions
collect `((nth ,i (storage-size object)) size)))
(class-name (data-class-name category cl-elt-type)))
`(progn
(defclass ,class-name (,(data-superclass-name category cl-elt-type))
((cl-elt-type :initform ',cl-elt-type :reader cl-elt-type
:allocation :class)))
(data-go ,class-name
,(or (member category '(matrix combination))))
(defmfun alloc ((object ,class-name))
,(gsl-name "alloc") ,cargs
:type :method
:c-return (cr :pointer)
:return ((assign-pointer object cr)))
(defmfun calloc ((object ,class-name))
,(gsl-name "calloc") ,cargs
:type :method
:c-return (cr :pointer)
:return ((assign-pointer object cr)))
(defmfun free ((object ,class-name))
,(gsl-name "free") (((pointer object) :pointer))
:type :method
:c-return :void)
(defmfun write-binary ((object ,class-name) stream)
,(gsl-name "fwrite") ((stream :pointer) ((pointer object) :pointer))
:type :method)
(defmfun read-binary ((object ,class-name) stream)
,(gsl-name "fread") ((stream :pointer) ((pointer object) :pointer))
:type :method)
(defmfun write-formatted ((object ,class-name) stream format)
,(gsl-name "fprintf")
((stream :pointer) ((pointer object) :pointer) (format :string))
:type :method)
(defmfun read-formatted ((object ,class-name) stream format)
,(gsl-name "fscanf")
((stream :pointer) ((pointer object) :pointer) (format :string))
:type :method)))))
(defun defmfun-all (category cl-types args &optional key-string)
"A defmfun for each of the declared data types."
(let ((categories (if (listp category) category (list category))))
`(progn
,@(loop for type in cl-types
for ctype = (lookup-C-type type)
collect
`(defmfun ,(first args)
;; Set the class name for the arglist
,(mapcar
(lambda (x)
(if (and (listp x) (member (second x) categories))
(list (first x)
(data-class-name (second x) type))
x))
(second args))
;; Create the correct GSL C library function name
,(splice-name
(third args)
;; Use key-string to override the lookup of the category.
(or key-string (data-gsl-string (first categories) type))
type)
,(mapcar
(lambda (x)
(if (eq (st-type x) :c-base-type)
(list (st-symbol x) ctype)
x))
(fourth args))
,@(let ((restargs (copy-list (nthcdr 4 args))))
(when (eq (getf restargs :c-return) :c-base-type)
(setf (getf restargs :c-return) ctype))
(setf (getf restargs :type) :method)
restargs))))))
;;;;****************************************************************************
;;;; Making data objects and initializing storage
;;;;****************************************************************************
(export '(make-data))
(defgeneric alloc (object)
(:documentation "Allocate GSL data; used internally."))
(defgeneric calloc (object)
(:documentation "Allocate GSL data and clear; used internally."))
(export 'free)
(defgeneric free (object)
(:documentation "Free GSL data.")
(:method :after ((object gsl-data)) (setf (pointer object) nil)))
(defun make-data (type zero &rest size)
"Make the GSL data object, including the allocation of space.
The user is responsible for calling #'free to free the foreign
memory when done."
(let ((obj (make-instance type :storage-size size)))
(if zero (calloc obj) (alloc obj))
obj))
;;;;****************************************************************************
;;;; Getting values into CL
;;;;****************************************************************************
(defun cl-invalidate (&rest objects)
"Mark the CL image of the GSL array as invalid."
(mapc (lambda (obj) (setf (cl-invalid obj) t))
objects))
(export 'data)
(defgeneric data (object &optional destination)
(:documentation "Extract the values in the object to a CL object.
The destination may be a sequence, 'vector, 'list.
If it is a sequence, that object is filled with the values.
If it is any other object, a new sequence is made of the type
specified.")
;; Default method is to make a sequence
(:method ((object gsl-data) &optional (sequence 'vector))
(if (eq (cl-invalid object) t)
;; set everything
(let* ((total-size (apply #'* (storage-size object)))
(seq
(case sequence
(list (make-list total-size))
((nil vector)
(make-array (list total-size)
:element-type (cl-elt-type object)))
(t sequence))))
(loop for i from 0
below (min (length seq) total-size)
do (setf (elt seq i) (maref object i)))
seq)
;; set selected
(let ((seq (data-cache object)))
(mapc (lambda (is)
(let ((i (first is)))
(setf (elt seq i) (maref object i))))
(cl-invalid object))
seq)))
;; Around method looks for cached value and returns it if valid;
;; otherwise computes CL element(s).
(:method :around ((object gsl-data) &optional (sequence 'vector))
(declare (ignore sequence))
(when (cl-invalid object)
(setf (data-cache object)
(call-next-method)
(cl-invalid object)
nil))
(data-cache object)))
;; Using GSL bulk data (vectors, matrices, etc.) storage.
;; Liam Healy, Sun Mar 26 2006 - 16:32
;; Time-stamp: <2008-04-26 22:13:16EDT mathematical.lisp>
;; $Id$
(in-package :gsl)
;;;;****************************************************************************
;;;; Setting values from CL
;;;;****************************************************************************
(defgeneric (setf data) (cl-array object)
(:documentation "Set the values in the object from a CL array.")
;; Default method is to read from a sequence
(:method (sequence (object gsl-data))
(loop for i from 0
below (min (length sequence) (apply #'* (storage-size object)))
do (setf (maref object i) (elt sequence i))))
(:method :after (source (object gsl-data))
(setf (cl-invalid object) t)))
(export 'set-identity)
(defgeneric set-identity (object)
(:documentation "Set elements to represent the identity.")
(:method :after ((object gsl-data)) (cl-invalidate object)))
(export 'set-basis)
(defgeneric set-basis (object index)
(:documentation "Set indexth basis vector."))
(export 'data-valid)
(defgeneric data-valid (object)
(:documentation "Validate the values in the object."))
;;;;****************************************************************************
;;;; Copying
;;;;****************************************************************************
(export '(copy swap))
(defgeneric copy (destination source)
(:documentation "Copy from source to destination."))
(defgeneric swap (obj1 obj2)
(:documentation "Swap contents of obj1 and obj2."))
;;;;****************************************************************************
;;;; Arithmetic operations
;;;;****************************************************************************
(export '(m+ m- m* m/ m*c m+c))
(defgeneric m+ (a b)
(:documentation "Add."))
(defgeneric m- (a b)
(:documentation "Subtract."))
(defgeneric m* (a b)
(:documentation "Multiply."))
(defgeneric m/ (a b)
(:documentation "Divide."))
(defgeneric m+c (a x)
(:documentation "Add scalar."))
(defgeneric m*c (a x)
(:documentation "Multiply by scalar."))
;;;;****************************************************************************
;;;; Maximum and minimum elements
;;;;****************************************************************************
(export
'(gsl-max gsl-min gsl-minmax gsl-max-index gsl-min-index gsl-minmax-index))
(defgeneric gsl-min (a)
(:documentation "Minimum."))
(defgeneric gsl-max (a)
(:documentation "Maximum."))
(defgeneric gsl-minmax (a)
(:documentation "Minimum and maximum."))
(defgeneric gsl-min-index (a)
(:documentation "Index of minimum."))
(defgeneric gsl-max-index (a)
(:documentation "Index of maximum."))
(defgeneric gsl-minmax-index (a)
(:documentation "Indices of minimum and maximum."))
;;;;****************************************************************************
;;;; Properties
;;;;****************************************************************************
(export '(gsl-zerop))
(defgeneric gsl-zerop (a)
(:documentation "Object is zero."))
;; Matrices
;; Liam Healy, Sun Mar 26 2006 - 11:51
;; Time-stamp: <2008-03-18 21:16:41EDT matrix.lisp>
;; $Id$
(in-package :gsl)
;;; Matrices are specified in a letm binding with
;;; (matrix-double-float size-or-initial &optional zero)
;;; (matrix-single-float size-or-initial &optional zero)
;;; (matrix-fixnum size-or-initial &optional zero)
;;; (matrix-complex size-or-initial &optional zero)
;;; where size-or-initial is a length-2 list of positive integers
;;; indicating the dimensions, and zero indicates that all elements should
;;; be set to zero, or, size-or-initial is a 2D array to which the
;;; vector should be initially set.
;;;;****************************************************************************
;;;; Matrix object definition, allocation, reading & writing
;;;;****************************************************************************
;;; GSL-matrix definition
(cffi:defcstruct gsl-matrix-c
(size1 size)
(size2 size)
(tda size)
(data :pointer)
(block :pointer)
(owner :int))
(defclass matrix (gsl-data) ())
(add-data-class matrix double-float matrix-double-float matrix "matrix")
(add-data-class matrix single-float matrix-single-float matrix "matrix")
(add-data-class matrix fixnum matrix-fixnum matrix "matrix")
(add-data-class matrix complex matrix-complex matrix "matrix")
(defdata matrix double-float 2)
(defdata matrix single-float 2)
(defdata matrix fixnum 2)
(defdata matrix complex 2)
(defmacro defmfun-mdsfc (&rest args)
"A defmfun for matrices of double, single, fixnum, and complex."
(defmfun-all '(matrix vector) '(double-float single-float fixnum complex) args))
(defmacro defmfun-mdsf (&rest args)
"A defmfun for matrices of double, single, and fixnum."
(defmfun-all '(matrix vector) '(double-float single-float fixnum) args))
(defmethod gsl-array ((object matrix))
(foreign-slot-value (pointer object) 'gsl-matrix-c 'data))
(export 'matrix-data)
(defun matrix-data (pointer)
"A pointer to the GSL array with the data contents, from the
sruct pointer."
(cffi:foreign-slot-value pointer 'gsl-matrix-c 'data))
;;;;****************************************************************************
;;;; Getting values
;;;;****************************************************************************
(defmfun-mdsfc maref ((matrix matrix) &rest indices)
"gsl_matrix_get"
(((pointer matrix) :pointer)
((first indices) size)
((second indices) size))
:c-return :c-base-type
:documentation ; FDL
"The (i,j)-th element of the matrix.")
(defmfun mref (pointer index0 index1)
"gsl_matrix_get"
((pointer :pointer) (index0 size) (index1 size))
:c-return :double
:index nil
:documentation ; FDL
"An element of the matrix of doubles, computed from the pointer.")
(export 'gsl-matrix-ptr)
(defgeneric gsl-matrix-ptr (matrix i j)
(:documentation
"A pointer to the i,j-th element of a matrix."))
(defmfun-mdsfc gsl-matrix-ptr ((matrix matrix) i j)
"gsl_matrix_ptr" (((pointer matrix) :pointer) (i size) (j size))
:c-return :pointer)
(defmethod data ((object matrix) &optional array)
(let ((arr (or array
(make-array (storage-size object)
:element-type (cl-elt-type object)))))
(loop for i from 0
below (min (array-dimension arr 0) (first (storage-size object)))
do
(loop for j from 0
below
(min (array-dimension arr 1) (second (storage-size object)))
do
(setf (aref arr i j) (maref object i j))))
arr))
;;;;****************************************************************************
;;;; Setting values
;;;;****************************************************************************
(defmfun-mdsfc (setf maref) (value (matrix matrix) &rest indices)
"gsl_matrix_set"
(((pointer matrix) :pointer)
((first indices) size)
((second indices) size)
(value :c-base-type))
:c-return :void
:documentation ; FDL
"Set the (i,j)-th element of the matrix.")
(defmfun (setf mref) (value pointer index0 index1)
"gsl_matrix_set"
((pointer :pointer) (index0 size) (index1 size) (value :double))
:c-return :void
:index nil
:documentation ; FDL
"Set an element of the matrix of doubles, using its pointer.")
(defmethod (setf data) (array (object matrix))
(loop for i from 0
below
(min (array-dimension array 0) (first (storage-size object)))
do
(loop for j from 0
below
(min (array-dimension array 1) (second (storage-size object)))
do
(setf (maref object i j) (aref array i j)))))
(defmfun-mdsfc set-all ((object matrix) value)
"gsl_matrix_set_all"
(((pointer object) :pointer) (value :c-base-type))
:c-return :void)
(defmfun-mdsfc set-zero ((object matrix))
"gsl_matrix_set_zero"
(((pointer object) :pointer))
:c-return :void)
(defmfun-mdsfc set-identity ((matrix matrix))
"gsl_matrix_set_identity" (((pointer matrix) gsl-matrix-c))
:c-return :void
:documentation ; FDL
"Set the elements of the matrix m to the
corresponding elements of the identity matrix, m(i,j) =
\delta(i,j), i.e. a unit diagonal with all off-diagonal elements zero.
This applies to both square and rectangular matrices.")
;;;;****************************************************************************
;;;; Matrix Views
;;;;****************************************************************************
(cffi:defcstruct gsl-matrix-view
(matrix gsl-matrix-c))
(export 'submatrix)
(defgeneric submatrix (matrix k1 k2 n1 n2)
(:documentation ; FDL
"A matrix view of a submatrix of the matrix.
The upper-left element of the submatrix is the element
(k1, k2) of the original matrix. The submatrix has n1
rows and n2 columns. The physical number of columns in memory
is unchanged."))
(defmfun-mdsfc submatrix ((matrix matrix) k1 k2 n1 n2)
"gsl_matrix_submatrix"
(((pointer matrix) gsl-matrix-c) (k1 size) (k2 size) (n1 size) (n2 size))
:c-return gsl-matrix-view)
(export 'matrix-array)
(defgeneric matrix-array (matrix n1 n2)
(:documentation ; FDL
"A matrix view of the array. The
matrix has n1 rows and n2 columns. The physical number of
columns in memory is also given by n2."))
(defmfun-mdsfc matrix-array ((matrix matrix) n1 n2)
"gsl_matrix_view_array"
(((pointer matrix) gsl-matrix-c) (n1 size) (n2 size))
:c-return gsl-matrix-view)
(export 'matrix-array-tda)
(defgeneric matrix-array-tda (matrix i j tda)
(:documentation ; FDL
"A matrix view of the array with a
physical number of columns tda which may differ from the corresponding
dimension of the matrix. The matrix has n1 rows and n2
columns, and the physical number of columns in memory is given by
tda."))
(defmfun-mdsfc matrix-array-tda ((matrix matrix) n1 n2 tda)
"gsl_matrix_view_array_with_tda"
(((pointer matrix) gsl-matrix-c) (n1 size) (n2 size) (tda size))
:c-return gsl-matrix-view)
(defmacro defmfun-mvdsfc (&rest args)
"A defmfun for vectors of double, single, fixnum, and complex,
translating to a GSL function named matrix_*."
(defmfun-all
'vector '(double-float single-float fixnum complex) args
"matrix"))
(export 'matrix-vector)
(defgeneric matrix-vector (vector n1 n2)
(:documentation ; FDL
"A matrix view of the vector. The matrix
has n1 rows and n2 columns. The vector must have unit
stride. The physical number of columns in memory is also given by
n2. Mathematically, the (i,j)-th element of the new
matrix is given by m'(i,j) = v->data[i*n2 + j]
where the index i runs from 0 to n1-1 and the index j
runs from 0 to n2-1.
The new matrix is only a view of the vector. When the view
goes out of scope the original vector will continue to exist.
The original memory can only be deallocated by freeing the original
vector. Of course, the original vector should not be deallocated while
the view is still in use."))
(defmfun-mvdsfc matrix-vector ((v vector) n1 n2)
"gsl_matrix_view_vector"
(((pointer v) gsl-vector-c) (n1 size) (n2 size))
:c-return gsl-matrix-view)
(export 'matrix-vector-tda)
(defgeneric matrix-vector-tda (vector n1 n2 tda)
(:documentation ; FDL
"A matrix view of the vector with a
physical number of columns tda which may differ from the
corresponding matrix dimension. The vector must have unit stride. The
matrix has n1 rows and n2 columns, and the physical number
of columns in memory is given by tda. Mathematically, the
(i,j)-th element of the new matrix is given by
m'(i,j) = v->data[i*tda + j]
where the index i runs from 0 to n1-1 and the index j
runs from 0 to n2-1.
The new matrix is only a view of the vector. When the view
goes out of scope the original vector will continue to exist.
The original memory can only be deallocated by freeing the original
vector. Of course, the original vector should not be deallocated while
the view is still in use."))
(defmfun-mvdsfc matrix-vector-tda ((v vector) n1 n2 tda)
"gsl_matrix_view_vector_with_tda"
(((pointer v) gsl-vector-c) (n1 size) (n2 size) (tda size))
:c-return gsl-matrix-view)
;;;;****************************************************************************
;;;; Row and Column Views
;;;;****************************************************************************
(export 'row-view)
(defgeneric row-view (matrix i)
(:documentation ; FDL
"A vector view of the ith row of the matrix."))
(defmfun-mdsfc row-view ((matrix matrix) i)
"gsl_matrix_row" (((pointer matrix) gsl-matrix-c) (i size))
:c-return gsl-vector-view
:null-pointer-info (:EFAULT (format nil "index ~d out of range" i)))
(export 'column-view)
(defgeneric column-view (matrix j)
(:documentation ; FDL
"A vector view of the jth column of the matrix."))
(defmfun-mdsfc column-view ((matrix matrix) j)
"gsl_matrix_column" ((matrix gsl-matrix-c) (j size))
:c-return gsl-matrix-view
:null-pointer-info (:EFAULT (format nil "index ~d out of range" j)))
(export 'diagonal-view)
(defgeneric diagonal-view (matrix)
(:documentation ; FDL
"A vector view of the diagonal of the matrix.
The matrix is not required to be square.
For a rectangular matrix the length of the diagonal is the same as the smaller
dimension of the matrix."))
(defmfun-mdsfc diagonal-view ((matrix matrix))
"gsl_matrix_diagonal" ((matrix gsl-matrix-c))
:c-return gsl-matrix-view)
(export 'subdiagonal-view)
(defgeneric subdiagonal-view (matrix k)
(:documentation ; FDL
"A vector view of the kth subdiagonal of the matrix; it is not
required to be square. The diagonal of the matrix corresponds to
k = 0."))
(defmfun-mdsfc subdiagonal-view ((matrix matrix) k)
"gsl_matrix_subdiagonal" ((matrix gsl-matrix-c) (k size))
:c-return gsl-matrix-view)
(export 'superdiagonal-view)
(defgeneric superdiagonal-view (matrix k)
(:documentation ; FDL
"A vector view of the kth superdiagonal of the matrix; it is not
required to be square. The diagonal of the matrix corresponds to k = 0."))
(defmfun-mdsfc superdiagonal-view ((matrix matrix) k)
"gsl_matrix_superdiagonal"
((matrix gsl-matrix-c) (k size))
:c-return gsl-matrix-view)
;;;;****************************************************************************
;;;; Copying
;;;;****************************************************************************
(defmfun-mdsfc copy ((destination matrix) (source matrix))
"gsl_matrix_memcpy"
(((pointer destination) gsl-matrix-c) ((pointer source) gsl-matrix-c))
:invalidate (destination)
:documentation ; FDL
"Copy the elements of the matrix source into the
matrix destination. The two matrices must have the same size.")
(defmfun-mdsfc swap ((m1 matrix) (m2 matrix))
"gsl_matrix_swap"
(((pointer m1) gsl-matrix-c) ((pointer m2) gsl-matrix-c))
:invalidate (m1 m2)
:documentation ; FDL
"Exchange the elements of the matrices m1 and
m2 by copying. The two matrices must have the same size.")
;;;;****************************************************************************
;;;; Copying rows and columns
;;;;****************************************************************************
(export 'row)
(defgeneric row (vector matrix i)
(:documentation ; FDL
"Copy the elements of the ith row of the matrix
into the vector. The length of the vector must be the
same as the length of the row."))
(defmfun-mdsfc row ((vector vector) (matrix matrix) i)
"gsl_matrix_get_row"
(((pointer vector) gsl-vector-c) ((pointer matrix) gsl-matrix-c) (i size))
:invalidate (vector))
(export 'column)
(defgeneric column (vector matrix j)
(:documentation ; FDL
"Copy the elements of the jth column of the matrix
into the vector. The length of the vector must be the
same as the length of the column."))
(defmfun-mdsfc column ((vector vector) (matrix matrix) j)
"gsl_matrix_get_col"
(((pointer vector) gsl-vector-c) ((pointer matrix) gsl-matrix-c) (j size))
:invalidate (vector))
(export 'set-row)
(defgeneric set-row (vector matrix i)
(:documentation ; FDL
"Copy the elements of the vector into the
ith row of the matrix. The length of the vector must be
the same as the length of the row."))
(defmfun-mdsfc set-row ((matrix matrix) i (vector vector))
"gsl_matrix_set_row"
(((pointer matrix) gsl-matrix-c) (i size) ((pointer vector) gsl-vector-c))
:invalidate (matrix))
(export 'set-column)
(defgeneric set-column (vector matrix j)
(:documentation ; FDL
"Copy the elements of the vector into the jth column of the matrix.
The length of the vector must be the same as the length of the column."))
(defmfun-mdsfc set-column ((matrix matrix) j (vector vector))
"gsl_matrix_set_col"
(((pointer matrix) gsl-matrix-c) (j size) ((pointer vector) gsl-vector-c))
:invalidate (matrix))
(defun (setf row) (vector matrix i) (set-row matrix i vector))
(defun (setf column) (vector matrix j) (set-column matrix j vector))
;;;;****************************************************************************
;;;; Exchanging rows and columns
;;;;****************************************************************************
(export 'swap-rows)
(defgeneric swap-rows (matrix i j)
(:documentation ; FDL
"Exchange the ith and jth rows of the matrix in-place."))
(defmfun-mdsfc swap-rows ((matrix matrix) i j)
"gsl_matrix_swap_rows"
(((pointer matrix) gsl-matrix-c) (i size) (j size))
:invalidate (matrix))
(export 'swap-columns)
(defgeneric swap-columns (matrix i j)
(:documentation ; FDL
"Exchange the ith and jth columns of the matrix in-place."))
(defmfun-mdsfc swap-columns ((matrix matrix) i j)
"gsl_matrix_swap_columns"
(((pointer matrix) gsl-matrix-c) (i size) (j size))
:invalidate (matrix))
(export 'swap-rowcol)
(defgeneric swap-rowcol (matrix i j)
(:documentation ; FDL
"Exchange the ith row and jth column of the
matrix in-place. The matrix must be square for this operation to
be possible."))
(defmfun-mdsfc swap-rowcol ((matrix matrix) i j)
"gsl_matrix_swap_rowcol"
(((pointer matrix) gsl-matrix-c) (i size) (j size))
:invalidate (matrix))
(export 'matrix-transpose-copy)
(defgeneric matrix-transpose-copy (destination source)
(:documentation ; FDL
"Make the destination matrix the transpose of the source matrix
by copying the elements. The dimensions of the destination
matrix must match the transposed dimensions of the source."))
(defmfun-mdsfc matrix-transpose-copy
((destination matrix) (source matrix))
"gsl_matrix_transpose_memcpy"
(((pointer destination) gsl-matrix-c) ((pointer source) gsl-matrix-c))
:invalidate (destination))
(export 'matrix-transpose)
(defgeneric matrix-transpose (matrix)
(:documentation ; FDL
"Replace the matrix by its transpose by copying the elements
of the matrix in-place. The matrix must be square for this
operation to be possible."))
(defmfun-mdsfc matrix-transpose ((matrix matrix))
"gsl_matrix_transpose"
(((pointer matrix) gsl-matrix-c))
:invalidate (matrix))
;;;;****************************************************************************
;;;; Arithmetic operations
;;;;****************************************************************************
(defmfun-mdsfc m+ ((a matrix) (b matrix))
"gsl_matrix_add"
(((pointer a) gsl-matrix-c) ((pointer b) gsl-matrix-c))
:invalidate (a)
:documentation ; FDL
"Add the elements of b to the elements of a,
a'_i = a_i + b_i. The two matrices must have the
same dimensions.")
(defmfun-mdsfc m- ((a matrix) (b matrix))
"gsl_matrix_sub" (((pointer a) gsl-matrix-c) ((pointer b) gsl-matrix-c))
:invalidate (a)
:documentation ; FDL
"Subtract the elements of matrix b from the elements of matrix
a, a'_i = a_i - b_i. The two matrices must have the
same dimensions.")
(defmfun-mdsfc m* ((a matrix) (b matrix))
"gsl_matrix_mul_elements"
(((pointer a) gsl-matrix-c) ((pointer b) gsl-matrix-c))
:invalidate (a)
:documentation ; FDL
"Multiply the elements of matrix a by the elements of
matrix b, a'(i,j) = a(i,j) * b(i,j). The two matrices must have the
same dimensions.")
(defmfun-mdsfc m/ ((a matrix) (b matrix))
"gsl_matrix_div_elements"
(((pointer a) gsl-matrix-c) ((pointer b) gsl-matrix-c))
:invalidate (a)
:documentation ; FDL
"Divide the elements of matrix a by the elements of
matrix b, a'(i,j) = a(i,j) / b(i,j). The two matrices must have the
same dimensions.")
(defmfun-mdsfc m*c ((a matrix) x)
"gsl_matrix_scale" (((pointer a) gsl-matrix-c) (x :c-base-type))
:invalidate (a)
:documentation ; FDL
"Multiply the elements of matrix a by the constant
factor x, a'(i,j) = x a(i,j).")
(defmfun-mdsfc m+c ((a matrix) x)
"gsl_matrix_add_constant" (((pointer a) gsl-matrix-c) (x :c-base-type))
:invalidate (a)
:documentation ; FDL
"Add the constant value x to the elements of the
matrix a, a'(i,j) = a(i,j) + x.")
;;;;****************************************************************************
;;;; Maximum and minimum elements
;;;;****************************************************************************
(defmfun-mdsf gsl-max ((m matrix))
"gsl_matrix_max" (((pointer m) gsl-matrix-c))
:documentation ; FDL
"The maximum value in the matrix m."
:c-return :c-base-type)
(defmfun-mdsf gsl-min ((m matrix))
"gsl_matrix_min" (((pointer m) gsl-matrix-c))
:documentation ; FDL
"The minimum value in the matrix m."
:c-return :c-base-type)
(defmfun-mdsf gsl-minmax ((m matrix))
"gsl_matrix_minmax"
(((pointer m) gsl-matrix-c) (min :c-base-type) (max :c-base-type))
:documentation ; FDL
"The minimum and maximum values in the matrix m."
:c-return :void)
(defmfun-mdsf gsl-max-index ((m matrix))
"gsl_matrix_max_index"
(((pointer m) gsl-matrix-c) (imax size) (jmax size))
:documentation ; FDL
"The index of the maximum value in the matrix m
When there are several equal maximum elements then the lowest index is
returned."
:c-return :void
:return ((list (scref imax) (scref jmax))))
(defmfun-mdsf gsl-min-index ((m matrix))
"gsl_matrix_min_index"
(((pointer m) gsl-matrix-c) (imin size) (jmin size))
:documentation ; FDL
"The index of the minimum value in the matrix m
When there are several equal minimum elements then the
lowest index is returned."
:c-return :void
:return ((list (scref imin) (scref jmin))))
(defmfun-mdsf gsl-minmax-index ((m matrix))
"gsl_matrix_minmax_index"
(((pointer m) gsl-matrix-c)
(imin size) (jmin size) (imax size) (jmax size))
:documentation ; FDL
"The indices of the minimum and maximum values in the matrix m.
When there are several equal minimum elements then the lowest index is
returned."
:c-return :void
:return ((list (scref imin) (scref jmin))
(list (scref imax) (scref jmax))))
;;;;****************************************************************************
;;;; Properties
;;;;****************************************************************************
(defmfun-mdsfc gsl-zerop ((m matrix))
"gsl_matrix_isnull" (((pointer m) gsl-matrix-c))
:documentation ; FDL
"All elements of matrix m are zero."
:c-return :boolean)
;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************
#|
(make-tests matrix-fixnum
(letm ((intmat (matrix-fixnum 2 2))) ;(setf maref), maref
(setf (maref intmat 0 1) 77)
(maref intmat 0 1))
(letm ((intmat (matrix-fixnum 2 2))) ;(setf data)
(setf (data intmat) #2A((4 6) (8 2)))
(data intmat))
(letm ((intmat (matrix-fixnum 2 2))) ;set-zero
(set-zero intmat)
(data intmat))
(letm ((intmat (matrix-fixnum 2 2))) ;set-all
(set-all intmat 44)
(data intmat))
(letm ((intmat (matrix-fixnum 2 2))) ;set-identity
(set-identity intmat)
(data intmat))
(letm ((intmat (matrix-fixnum #2A((4 6) (8 2)))) ;row
(vect (vector-fixnum 2)))
(row vect intmat 0)
(data vect))
(letm ((intmat (matrix-fixnum #2A((4 6) (8 2)))) ;column
(vect (vector-fixnum 2)))
(column vect intmat 1)
(data vect))
(letm ((intmat (matrix-fixnum #2A((-1 -12) (8 3))))) ;gsl-min
(gsl-min intmat))
(letm ((intmat (matrix-fixnum #2A((-1 -12) (8 3))))) ;gsl-max
(gsl-max intmat))
(letm ((intmat (matrix-fixnum #2A((-1 -12) (8 3))))) ;gsl-minmax
(multiple-value-list (gsl-minmax intmat)))
(letm ((intmat (matrix-fixnum #2A((-1 -12) (8 3))))) ;gsl-min-index
(gsl-min-index intmat))
(letm ((intmat (matrix-fixnum #2A((-1 -12) (8 3))))) ;gsl-max-index
(gsl-max-index intmat))
(letm ((intmat (matrix-fixnum #2A((-1 -12) (8 3))))) ;gsl-minmax-index
(multiple-value-list (gsl-minmax-index intmat)))
(letm ((intmat1 (matrix-fixnum #2A((1 2)(3 4)))) ;copy
(intmat2 (matrix-fixnum 2 2)))
(copy intmat2 intmat1)
(data intmat2))
(letm ((intmat1 (matrix-fixnum #2A((1 2)(3 4)))) ;swap
(intmat2 (matrix-fixnum #2A((5 6) (7 8)))))
(swap intmat1 intmat2)
(data intmat1))
(letm ((intmat1 (matrix-fixnum #2A((1 2)(3 4))))) ;swap-rows
(swap-rows intmat1 0 1)
(data intmat1))
(letm ((intmat1 (matrix-fixnum #2A((1 2)(3 4))))) ;swap-columns
(swap-columns intmat1 0 1)
(data intmat1))
(letm ((intmat1 (matrix-fixnum #2A((1 2)(3 4))))) ;swap-rowcol
(swap-rowcol intmat1 0 1)
(data intmat1)))
|#
(LISP-UNIT:DEFINE-TEST
MATRIX-FIXNUM
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 77)
(MULTIPLE-VALUE-LIST
(LETM ((INTMAT (MATRIX-FIXNUM 2 2)))
(SETF (MAREF INTMAT 0 1) 77)
(MAREF INTMAT 0 1))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #2A((4 6) (8 2)))
(MULTIPLE-VALUE-LIST
(LETM ((INTMAT (MATRIX-FIXNUM 2 2)))
(SETF (DATA INTMAT) #2A((4 6) (8 2)))
(DATA INTMAT))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #2A((0 0) (0 0)))
(MULTIPLE-VALUE-LIST
(LETM ((INTMAT (MATRIX-FIXNUM 2 2))) (SET-ZERO INTMAT)
(DATA INTMAT))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #2A((44 44) (44 44)))
(MULTIPLE-VALUE-LIST
(LETM ((INTMAT (MATRIX-FIXNUM 2 2)))
(SET-ALL INTMAT 44) (DATA INTMAT))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #2A((1 0) (0 1)))
(MULTIPLE-VALUE-LIST
(LETM ((INTMAT (MATRIX-FIXNUM 2 2)))
(SET-IDENTITY INTMAT) (DATA INTMAT))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #(4 6))
(MULTIPLE-VALUE-LIST
(LETM ((INTMAT (MATRIX-FIXNUM #2A((4 6) (8 2))))
(VECT (VECTOR-FIXNUM 2)))
(ROW VECT INTMAT 0)
(DATA VECT))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #(6 2))
(MULTIPLE-VALUE-LIST
(LETM ((INTMAT (MATRIX-FIXNUM #2A((4 6) (8 2))))
(VECT (VECTOR-FIXNUM 2)))
(COLUMN VECT INTMAT 1)
(DATA VECT))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST -12)
(MULTIPLE-VALUE-LIST
(LETM ((INTMAT (MATRIX-FIXNUM #2A((-1 -12) (8 3)))))
(GSL-MIN INTMAT))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 8)
(MULTIPLE-VALUE-LIST
(LETM ((INTMAT (MATRIX-FIXNUM #2A((-1 -12) (8 3)))))
(GSL-MAX INTMAT))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST (LIST -12 8))
(MULTIPLE-VALUE-LIST
(LETM ((INTMAT (MATRIX-FIXNUM #2A((-1 -12) (8 3)))))
(MULTIPLE-VALUE-LIST (GSL-MINMAX INTMAT)))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST (LIST 0 1))
(MULTIPLE-VALUE-LIST
(LETM ((INTMAT (MATRIX-FIXNUM #2A((-1 -12) (8 3)))))
(GSL-MIN-INDEX INTMAT))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST (LIST 1 0))
(MULTIPLE-VALUE-LIST
(LETM ((INTMAT (MATRIX-FIXNUM #2A((-1 -12) (8 3)))))
(GSL-MAX-INDEX INTMAT))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST (LIST (LIST 0 1) (LIST 1 0)))
(MULTIPLE-VALUE-LIST
(LETM ((INTMAT (MATRIX-FIXNUM #2A((-1 -12) (8 3)))))
(MULTIPLE-VALUE-LIST
(GSL-MINMAX-INDEX INTMAT)))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #2A((1 2) (3 4)))
(MULTIPLE-VALUE-LIST
(LETM ((INTMAT1 (MATRIX-FIXNUM #2A((1 2) (3 4))))
(INTMAT2 (MATRIX-FIXNUM 2 2)))
(COPY INTMAT2 INTMAT1) (DATA INTMAT2))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #2A((5 6) (7 8)))
(MULTIPLE-VALUE-LIST
(LETM ((INTMAT1 (MATRIX-FIXNUM #2A((1 2) (3 4))))
(INTMAT2 (MATRIX-FIXNUM #2A((5 6) (7 8)))))
(SWAP INTMAT1 INTMAT2) (DATA INTMAT1))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #2A((3 4) (1 2)))
(MULTIPLE-VALUE-LIST
(LETM ((INTMAT1 (MATRIX-FIXNUM #2A((1 2) (3 4)))))
(SWAP-ROWS INTMAT1 0 1) (DATA INTMAT1))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #2A((2 1) (4 3)))
(MULTIPLE-VALUE-LIST
(LETM ((INTMAT1 (MATRIX-FIXNUM #2A((1 2) (3 4)))))
(SWAP-COLUMNS INTMAT1 0 1) (DATA INTMAT1))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #2A((2 4) (3 1)))
(MULTIPLE-VALUE-LIST
(LETM ((INTMAT1 (MATRIX-FIXNUM #2A((1 2) (3 4)))))
(SWAP-ROWCOL INTMAT1 0 1) (DATA INTMAT1)))))
#|
(make-tests
matrix-double
(letm ((mat (matrix-double-float 10 3)))
(loop for i from 0 below 10
do
(loop for j from 0 below 3
do (setf (maref mat i j) (+ 0.23d0 j (* 100 i)))))
(data mat))
(letm ((mat (matrix-double-float #2A((1.0d0 2.0d0) (3.0d0 4.0d0))))
(ans (matrix-double-float 2 2)))
(copy ans mat)
(data ans))
(letm ((mat (matrix-double-float #2A((1.0d0 2.0d0) (3.0d0 4.0d0)))))
(m* mat mat)
(data mat))
(letm ((mat (matrix-double-float #2A((1.0d0 2.0d0) (-3.0d0 4.0d0)))))
(list (maref mat 0 0) (maref mat 0 1) (maref mat 1 0) (maref mat 1 1))))
|#
(LISP-UNIT:DEFINE-TEST MATRIX-DOUBLE
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST
#2A((0.23d0 1.23d0 2.23d0)
(100.23d0 101.23d0 102.23d0)
(200.23d0 201.23d0 202.23d0)
(300.23d0 301.23d0 302.23d0)
(400.23d0 401.23d0 402.23d0)
(500.23d0 501.23d0 502.23d0)
(600.23d0 601.23d0 602.23d0)
(700.23d0 701.23d0 702.23d0)
(800.23d0 801.23d0 802.23d0)
(900.23d0 901.23d0 902.23d0)))
(MULTIPLE-VALUE-LIST
(LETM ((MAT (MATRIX-DOUBLE-FLOAT 10 3)))
(LOOP FOR I FROM 0 BELOW 10 DO
(LOOP FOR J FROM 0 BELOW 3 DO
(SETF (MAREF MAT I J)
(+ 0.23d0 J (* 100 I)))))
(DATA MAT))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #2A((1.0d0 2.0d0) (3.0d0 4.0d0)))
(MULTIPLE-VALUE-LIST
(LETM ((MAT (MATRIX-DOUBLE-FLOAT #2A((1.0d0 2.0d0) (3.0d0 4.0d0))))
(ANS (MATRIX-DOUBLE-FLOAT 2 2)))
(COPY ANS MAT) (DATA ANS))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #2A((1.0d0 4.0d0) (9.0d0 16.0d0)))
(MULTIPLE-VALUE-LIST
(LETM ((MAT (MATRIX-DOUBLE-FLOAT #2A((1.0d0 2.0d0) (3.0d0 4.0d0)))))
(M* MAT MAT) (DATA MAT))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST (LIST 1.0d0 2.0d0 -3.0d0 4.0d0))
(MULTIPLE-VALUE-LIST
(LETM ((MAT (MATRIX-DOUBLE-FLOAT #2A((1.0d0 2.0d0) (-3.0d0 4.0d0)))))
(LIST (MAREF MAT 0 0) (MAREF MAT 0 1) (MAREF MAT 1 0) (MAREF MAT 1 1))))))
;; Vectors
;; Liam Healy, Sun Mar 26 2006 - 11:51
;; Time-stamp: <2008-04-13 09:38:52EDT vector.lisp>
;; $Id$
(in-package :gsl)
;;; Vectors are specified in a letm binding with
;;; (vector-double-float size-or-initial &optional zero)
;;; (vector-single-float size-or-initial &optional zero)
;;; (vector-fixnum size-or-initial &optional zero)
;;; (vector-complex size-or-initial &optional zero)
;;; where size-or-initial is a positive integer indicating the
;;; size, and zero indicates that all elements should be set to zero,
;;; or, size-or-initial is a sequence to which the vector
;;; should be initially set.
;;; No mechanism for C stream input/output yet.
;;; Generalize check-gsl-status to optionally signal errors, use here?
;;; Functions like write-binary etc. as a single function, selecting the C fn with typecase?
;;; #'subvector, #'subvector-stride cause crash, see notes 2006-03-30
;;; #'vector-complex-real, #'vector-complex-imag need structure definition
;;; GSL bug?: no gsl_vector_complex_add, etc. ===> no "vector"
;;; Need to build real vector out of view pointer.
;;; Need to #'cl-invalidate when (setf maref) called, see Mon Nov 26 2007.
;;;;****************************************************************************
;;;; Vector structure and CL object
;;;;****************************************************************************
#|
;;FDL
The size is simply the number of vector elements. The range of
valid indices runs from 0 to size-1. The stride is the
step-size from one element to the next in physical memory, measured in
units of the appropriate datatype. The pointer data gives the
location of the first element of the vector in memory. The pointer
block stores the location of the memory block in which the vector
elements are located (if any). If the vector owns this block then the
owner field is set to one and the block will be deallocated when the
vector is freed. If the vector points to a block owned by another
object then the owner field is zero and any underlying block will not be
deallocated with the vector.
|#
;;; GSL-vector definition
(cffi:defcstruct gsl-vector-c
(size size)
(stride size)
(data :pointer)
(block :pointer)
(owner :int))
(defclass mvector (gsl-data) ())
(add-data-class vector double-float vector-double-float mvector "vector")
(add-data-class vector single-float vector-single-float mvector "vector")
(add-data-class vector fixnum vector-fixnum mvector "vector")
(add-data-class vector complex vector-complex mvector "vector")
;;; Allocation, freeing, reading and writing
(defdata vector double-float)
(defdata vector single-float)
(defdata vector fixnum)
(defdata vector complex)
(defmacro defmfun-vdsfc (&rest args)
"A defmfun for vectors of double, single, fixnum, complex, unsigned-fixnum."
(defmfun-all 'vector '(double-float single-float fixnum complex) args))
(defmacro defmfun-vdsf (&rest args)
"A defmfun for vectors of double, single, and fixnum."
(defmfun-all 'vector '(double-float single-float fixnum) args))
(defmethod gsl-array ((object mvector))
(cffi:foreign-slot-value (pointer object) 'gsl-vector-c 'data))
(defun make-data-from-pointer (pointer &optional (class 'vector-double-float) size)
"Given a C pointer to a GSL data type, make the CL object."
(make-instance
class
:pointer pointer
:storage-size
(or size (cffi:foreign-slot-value pointer 'gsl-vector-c 'size))))
(export 'vector-data)
(defun vector-data (pointer)
"A pointer to the GSL array with the data contents, from the
sruct pointer."
(cffi:foreign-slot-value pointer 'gsl-vector-c 'data))
(defun element-size (object)
"The size of each element as stored in C."
(cffi:foreign-type-size
(rest (assoc (cl-elt-type object)
'((double . :double) (single . :float)
(fixnum . :int) (complex . gsl-complex))))))
;;;;****************************************************************************
;;;; Getting values
;;;;****************************************************************************
(defmfun-vdsfc maref ((vector vector) &rest indices)
"gsl_vector_get"
(((pointer vector) :pointer) ((first indices) size))
:c-return :c-base-type
:documentation ; FDL
"The ith element of the vector.")
(defmfun vref (pointer index)
"gsl_vector_get"
((pointer :pointer) (index size))
:c-return :double
:index nil
:documentation "An element of the vector of doubles, computed from the pointer.")
(defmfun gsl-vector-ptr (vector i)
"gsl_vector_ptr" (((pointer vector) :pointer) (i size))
:c-return :pointer
:documentation ; FDL
"The ith element of the vector as a pointer.")
;;;;****************************************************************************
;;;; Setting values
;;;;****************************************************************************
(defmfun-vdsfc (setf maref) (value (vector vector) &rest indices)
"gsl_vector_set"
(((pointer vector) :pointer) ((first indices) size) (value :c-base-type))
:c-return :void
:documentation ; FDL
"Set an element of the vector.")
(defmfun (setf vref) (value pointer index)
"gsl_vector_set"
((pointer :pointer) (index size) (value :double))
:c-return :void
:index nil
:documentation ; FDL
"Set an element of the vector of doubles, using its pointer.")
(defmfun-vdsfc set-all ((object vector) value)
"gsl_vector_set_all"
(((pointer object) :pointer) (value :c-base-type))
:c-return :void)
(defmfun-vdsfc set-zero ((object vector))
"gsl_vector_set_zero" (((pointer object) :pointer))
:c-return :void)
(defmfun-vdsfc set-basis ((vector vector) index)
"gsl_vector_set_basis" (((pointer vector) gsl-vector-c) (index size))
:invalidate (vector)
:documentation ; FDL
"Set the index element to 1, and the rest to 0.")
;;;;****************************************************************************
;;;; Views
;;;;****************************************************************************
(cffi:defcstruct gsl-vector-view
(vector gsl-vector-c))
;;; broken
;;; (letm ((vec (vector-double-float #(-3.21d0 1.0d0 12.8d0)))) (subvector vec 1 2))
(export '(subvector subvector-stride))
(defgeneric subvector (vector offset size)
(:documentation ; FDL
"Return a vector view of a subvector of another vector
v. The start of the new vector is offset by offset elements
from the start of the original vector. The new vector has size
elements."))
(defgeneric subvector-stride (vector offset stride size)
(:documentation ; FDL
"A vector view of a subvector of another vector
v with an additional stride argument. The subvector is formed in
the same way as for #'subvector but the new vector has
n elements with a step-size of stride from one element to
the next in the original vector. Mathematically, the i-th element
of the new vector v' is given by
v'(i) = v->data[(offset + i*stride)*v->stride]
where the index i runs from 0 to n-1.
Note that subvector views give direct access to the underlying elements
of the original vector."))
(defmfun-vdsfc subvector ((vector vector) offset size)
"gsl_vector_subvector"
(((pointer vector) gsl-vector-c) (offset size) (size size))
:c-return :pointer)
(defmfun-vdsfc subvector-stride ((vector vector) offset stride size)
"gsl_vector_subvector_with_stride"
(((pointer vector) gsl-vector-c)
(offset size) (stride size) (size size))
:c-return :pointer)
;;; These require that the gsl-vector-complex structure be defined.
#|
(defmfun vector-complex-real (vector)
"gsl_vector_complex_real" ((vector gsl-vector-complex))
:c-return :pointer
:documentation
"A vector view of the real parts of the complex vector v.")
(defmfun vector-complex-imag (vector)
"gsl_vector_complex_imag"((vector gsl-vector-complex))
:c-return :pointer
:documentation
"A vector view of the imaginary parts of the complex vector v.")
|#
(defmfun vector-array (base size)
"gsl_vector_view_array" ((base :pointer) (size size))
:c-return :pointer
:documentation ; FDL
"A vector view of an array. The start of the new
vector is given by base and has n elements.")
(defmfun vector-array-stride (base stride size)
"gsl_vector_view_array_with_stride"
((base :pointer) (stride size) (size size))
:c-return :pointer
:documentation ; FDL
"A vector view of an array with stride. The start of the new
vector is given by base.")
;;;;****************************************************************************
;;;; Copying
;;;;****************************************************************************
(defmfun-vdsfc copy ((destination vector) (source vector))
"gsl_vector_memcpy"
(((pointer destination) gsl-vector-c) ((pointer source) gsl-vector-c))
:invalidate (destination)
:documentation ; FDL
"Copy the elements of the vector source into the
vector destination. The two vectors must have the same length.")
(defmfun-vdsfc swap ((v vector) (w vector))
"gsl_vector_swap" (((pointer v) gsl-vector-c) ((pointer w) gsl-vector-c))
:invalidate (v w)
:documentation ; FDL
"Exchange the elements of the vectors v and w
by copying. The two vectors must have the same length.")
;;;;****************************************************************************
;;;; Exchanging elements
;;;;****************************************************************************
(export '(swap-elements vector-reverse))
(defgeneric swap-elements (vec i j)
(:documentation ; FDL
"Exchange the i-th and j-th elements of the vector vec in-place."))
(defgeneric vector-reverse (vec)
(:documentation
"Exchange the i-th and j-th elements of the vector vec in-place."))
(defmfun-vdsfc swap-elements ((vec vector) i j)
"gsl_vector_swap_elements" (((pointer vec) gsl-vector-c) (i size) (j size))
:after ((when (listp (cl-invalid vec))
(push (list i) (cl-invalid vec))
(push (list j) (cl-invalid vec))))
:return (vec))
(defmfun-vdsfc vector-reverse ((vec vector))
"gsl_vector_reverse" (((pointer vec) gsl-vector-c))
:invalidate (vec)
:documentation ; FDL
"Reverse the order of the elements of the vector vec.")
;;;;****************************************************************************
;;;; Arithmetic operations
;;;;****************************************************************************
(defmfun-vdsf m+ ((a vector) (b vector))
"gsl_vector_add" (((pointer a) gsl-vector-c) ((pointer b) gsl-vector-c))
:invalidate (a)
:documentation ; FDL
"Add the elements of vector b to the elements of
vector a, a'_i = a_i + b_i. The two vectors must have the
same length.")
(defmfun-vdsf m- ((a vector) (b vector))
"gsl_vector_sub" (((pointer a) gsl-vector-c) ((pointer b) gsl-vector-c))
:invalidate (a)
:documentation ; FDL
"Subtract the elements of vector b from the elements of
vector a, a'_i = a_i - b_i. The two vectors must have the same length.")
(defmfun-vdsf m* ((a vector) (b vector))
"gsl_vector_mul" (((pointer a) gsl-vector-c) ((pointer b) gsl-vector-c))
:invalidate (a)
:documentation ; FDL
"Multiply the elements of vector a by the elements of
vector b, a'_i = a_i * b_i. The two vectors must have the same length.")
(defmfun-vdsf m/ ((a vector) (b vector))
"gsl_vector_div" (((pointer a) gsl-vector-c) ((pointer b) gsl-vector-c))
:invalidate (a)
:documentation ; FDL
"Divide the elements of vector a by the elements of
vector b, a'_i = a_i / b_i. The two vectors must have the same length.")
(defmfun-vdsf m*c ((a vector) x)
"gsl_vector_scale" (((pointer a) gsl-vector-c) (x :double))
:invalidate (a)
:documentation ; FDL
"Multiply the elements of vector a by the constant factor x, a'_i = x a_i.")
(defmfun-vdsf m+c ((a vector) x)
"gsl_vector_add_constant" (((pointer a) gsl-vector-c) (x :double))
:invalidate (a)
:documentation ; FDL
"Add the constant value x to the elements of the vector a, a'_i = a_i + x.")
;;;;****************************************************************************
;;;; Maximum and minimum elements
;;;;****************************************************************************
(defmfun-vdsf gsl-max ((v vector))
"gsl_vector_max" (((pointer v) gsl-vector-c))
:c-return :c-base-type
:documentation ; FDL
"The maximum value in the vector v.")
(defmfun-vdsf gsl-min ((v vector))
"gsl_vector_min" (((pointer v) gsl-vector-c))
:c-return :c-base-type
:documentation ; FDL
"The minimum value in the vector v.")
(defmfun-vdsf gsl-minmax ((v vector))
"gsl_vector_minmax"
(((pointer v) gsl-vector-c) (min :c-base-type) (max :c-base-type))
:c-return :void
:documentation ; FDL
"The minimum and maximum values in the vector v.")
(defmfun-vdsf gsl-max-index ((v vector))
"gsl_vector_max_index" (((pointer v) gsl-vector-c))
:c-return size
:documentation ; FDL
"The index of the maximum value in the vector v.
When there are several equal minimum elements then the lowest index is
returned.")
(defmfun-vdsf gsl-min-index ((v vector))
"gsl_vector_min_index" (((pointer v) gsl-vector-c))
:c-return size
:documentation ; FDL
"The index of the minimum value in the vector v. When there are several
equal minimum elements then the lowest index is returned.")
(defmfun-vdsf gsl-minmax-index ((v vector))
"gsl_vector_minmax_index"
(((pointer v) gsl-vector-c) (imin size) (imax size))
:c-return :void
:documentation ; FDL
"The indices of the minimum and maximum values in the vector v.
When there are several equal minimum elements then the lowest index is
returned.")
;;;;****************************************************************************
;;;; Properties
;;;;****************************************************************************
(defmfun-vdsfc gsl-zerop ((v vector))
"gsl_vector_isnull" (((pointer v) gsl-vector-c))
:c-return :boolean
:documentation ; FDL
"All elements of vector v are zero.")
;;;;****************************************************************************
;;;; Examples and unit tests
;;;;****************************************************************************
#|
(make-tests
vector-fixnum
(letm ((intvec (vector-fixnum 4))) ;(setf maref), maref
(setf (maref intvec 1) 77)
(maref intvec 1))
(letm ((intvec (vector-fixnum 4))) ;(setf data)
(setf (data intvec) #(4 6 8 2))
(data intvec))
(letm ((intvec (vector-fixnum 4))) ;set-zero
(set-zero intvec)
(data intvec))
(letm ((intvec (vector-fixnum 4))) ;set-all
(set-all intvec 44)
(data intvec))
(letm ((intvec (vector-fixnum 4))) ;set-basis
(set-basis intvec 1)
(data intvec))
(letm ((intvec (vector-fixnum #(1 2 3 4)))) ;vector-reverse
(vector-reverse intvec)
(data intvec))
(letm ((intvec (vector-fixnum #(-1 -12 8 3)))) ;gsl-min
(gsl-min intvec))
(letm ((intvec (vector-fixnum #(-1 -12 8 3)))) ;gsl-max
(gsl-max intvec))
(letm ((intvec (vector-fixnum #(-1 -12 8 3)))) ;gsl-minmax
(multiple-value-list (gsl-minmax intvec)))
(letm ((intvec (vector-fixnum #(-1 -12 8 3)))) ;gsl-min-index
(gsl-min-index intvec))
(letm ((intvec (vector-fixnum #(-1 -12 8 3)))) ;gsl-max-index
(gsl-max-index intvec))
(letm ((intvec (vector-fixnum #(-1 -12 8 3)))) ;gsl-minmax-index
(multiple-value-list (gsl-minmax-index intvec)))
(letm ((intvec1 (vector-fixnum #(1 2 3 4))) ;copy
(intvec2 (vector-fixnum 4)))
(copy intvec2 intvec1)
(data intvec2))
(letm ((intvec1 (vector-fixnum #(1 2 3 4))) ;swap
(intvec2 (vector-fixnum #(5 6 7 8))))
(swap intvec2 intvec1)
(concatenate 'vector (data intvec1) (data intvec2)))
(letm ((intvec (vector-fixnum #(1 2 3 4)))) ;swap-elements
(swap-elements intvec 1 3)
(data intvec)))
|#
(LISP-UNIT:DEFINE-TEST VECTOR-FIXNUM
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 77)
(MULTIPLE-VALUE-LIST
(LETM ((INTVEC (VECTOR-FIXNUM 4)))
(SETF (MAREF INTVEC 1) 77)
(MAREF INTVEC 1))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #(4 6 8 2))
(MULTIPLE-VALUE-LIST
(LETM ((INTVEC (VECTOR-FIXNUM 4)))
(SETF (DATA INTVEC) #(4 6 8 2))
(DATA INTVEC))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #(0 0 0 0))
(MULTIPLE-VALUE-LIST
(LETM ((INTVEC (VECTOR-FIXNUM 4)))
(SET-ZERO INTVEC)
(DATA INTVEC))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #(44 44 44 44))
(MULTIPLE-VALUE-LIST
(LETM ((INTVEC (VECTOR-FIXNUM 4)))
(SET-ALL INTVEC 44)
(DATA INTVEC))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #(0 1 0 0))
(MULTIPLE-VALUE-LIST
(LETM ((INTVEC (VECTOR-FIXNUM 4)))
(SET-BASIS INTVEC 1)
(DATA INTVEC))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #(4 3 2 1))
(MULTIPLE-VALUE-LIST
(LETM ((INTVEC (VECTOR-FIXNUM #(1 2 3 4))))
(VECTOR-REVERSE INTVEC)
(DATA INTVEC))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST -12)
(MULTIPLE-VALUE-LIST
(LETM ((INTVEC (VECTOR-FIXNUM #(-1 -12 8 3))))
(GSL-MIN INTVEC))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 8)
(MULTIPLE-VALUE-LIST
(LETM ((INTVEC (VECTOR-FIXNUM #(-1 -12 8 3))))
(GSL-MAX INTVEC))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST (LIST -12 8))
(MULTIPLE-VALUE-LIST
(LETM ((INTVEC (VECTOR-FIXNUM #(-1 -12 8 3))))
(MULTIPLE-VALUE-LIST
(GSL-MINMAX INTVEC)))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 1)
(MULTIPLE-VALUE-LIST
(LETM ((INTVEC (VECTOR-FIXNUM #(-1 -12 8 3))))
(GSL-MIN-INDEX INTVEC))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 2)
(MULTIPLE-VALUE-LIST
(LETM
((INTVEC (VECTOR-FIXNUM #(-1 -12 8 3))))
(GSL-MAX-INDEX INTVEC))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST (LIST 1 2))
(MULTIPLE-VALUE-LIST
(LETM ((INTVEC (VECTOR-FIXNUM #(-1 -12 8 3))))
(MULTIPLE-VALUE-LIST
(GSL-MINMAX-INDEX INTVEC)))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #(1 2 3 4))
(MULTIPLE-VALUE-LIST
(LETM ((INTVEC1 (VECTOR-FIXNUM #(1 2 3 4)))
(INTVEC2 (VECTOR-FIXNUM 4)))
(COPY INTVEC2 INTVEC1)
(DATA INTVEC2))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #(5 6 7 8 1 2 3 4))
(MULTIPLE-VALUE-LIST
(LETM ((INTVEC1 (VECTOR-FIXNUM #(1 2 3 4)))
(INTVEC2 (VECTOR-FIXNUM #(5 6 7 8))))
(SWAP INTVEC2 INTVEC1)
(CONCATENATE 'VECTOR (DATA INTVEC1) (DATA INTVEC2)))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #(1 4 3 2))
(MULTIPLE-VALUE-LIST
(LETM ((INTVEC (VECTOR-FIXNUM #(1 2 3 4))))
(SWAP-ELEMENTS INTVEC 1 3)
(DATA INTVEC)))))
#|
(make-tests
vector-double
(letm ((vec (vector-double-float 3)))
(setf (maref vec 0) -3.21d0
(maref vec 1) 1.0d0
(maref vec 2) 12.8d0
(cl-invalid vec) t)
(data vec))
(letm ((vec (vector-double-float 3)))
(setf (data vec) #(-3.21d0 1.0d0 12.8d0))
(data vec))
(letm ((vec (vector-double-float #(-3.21d0 1.0d0 12.8d0))))
(data vec))
(letm ((base (vector-double-float 5)))
(set-basis base 1)
(data base))
(letm ((vec1 (vector-double-float #(-3.21d0 1.0d0 12.8d0)))
(vec2 (vector-double-float 3)))
(copy vec2 vec1)
(data vec2)))
|#
(LISP-UNIT:DEFINE-TEST VECTOR-DOUBLE-FLOAT
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #(-3.21d0 1.0d0 12.8d0))
(MULTIPLE-VALUE-LIST
(LETM ((VEC (VECTOR-DOUBLE-FLOAT 3)))
(SETF (MAREF VEC 0) -3.21d0
(MAREF VEC 1) 1.0d0
(MAREF VEC 2) 12.8d0
(CL-INVALID VEC) T)
(DATA VEC))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #(-3.21d0 1.0d0 12.8d0))
(MULTIPLE-VALUE-LIST
(LETM ((VEC (VECTOR-DOUBLE-FLOAT 3)))
(SETF (DATA VEC) #(-3.21d0 1.0d0 12.8d0))
(DATA VEC))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #(-3.21d0 1.0d0 12.8d0))
(MULTIPLE-VALUE-LIST
(LETM ((VEC (VECTOR-DOUBLE-FLOAT #(-3.21d0 1.0d0 12.8d0))))
(DATA VEC))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #(0.0d0 1.0d0 0.0d0 0.0d0 0.0d0))
(MULTIPLE-VALUE-LIST
(LETM ((BASE (VECTOR-DOUBLE-FLOAT 5))) (SET-BASIS BASE 1)
(DATA BASE))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST #(-3.21d0 1.0d0 12.8d0))
(MULTIPLE-VALUE-LIST
(LETM
((VEC1 (VECTOR-DOUBLE-FLOAT #(-3.21d0 1.0d0 12.8d0)))
(VEC2 (VECTOR-DOUBLE-FLOAT 3)))
(COPY VEC2 VEC1) (DATA VEC2)))))
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment