diff --git a/data/block.lisp b/data/block.lisp deleted file mode 100644 index 1d9a600a5853573136f9b186c03021e0e7138b35..0000000000000000000000000000000000000000 --- a/data/block.lisp +++ /dev/null @@ -1,21 +0,0 @@ -;; 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) diff --git a/data/data-expansions.lisp b/data/data-expansions.lisp deleted file mode 100644 index 23786ceabeddf81381194f8da88b7b7b66190dce..0000000000000000000000000000000000000000 --- a/data/data-expansions.lisp +++ /dev/null @@ -1,11 +0,0 @@ -;; 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." - ) diff --git a/data/data.lisp b/data/data.lisp deleted file mode 100644 index a73b911f352ed811d57629c9050a4df2adb29c92..0000000000000000000000000000000000000000 --- a/data/data.lisp +++ /dev/null @@ -1,335 +0,0 @@ -;; 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))) diff --git a/data/mathematical.lisp b/data/mathematical.lisp deleted file mode 100644 index df8c2d8cacf39877ee2e2e53ffa233ec535ea080..0000000000000000000000000000000000000000 --- a/data/mathematical.lisp +++ /dev/null @@ -1,102 +0,0 @@ -;; 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.")) diff --git a/data/matrix.lisp b/data/matrix.lisp deleted file mode 100644 index 669dfc089ab64a07ff2df663dec1807990cfe405..0000000000000000000000000000000000000000 --- a/data/matrix.lisp +++ /dev/null @@ -1,783 +0,0 @@ -;; 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)))))) - diff --git a/data/vector.lisp b/data/vector.lisp deleted file mode 100644 index c9de721f4abe708e41d2740b48ca26e195dbb4a6..0000000000000000000000000000000000000000 --- a/data/vector.lisp +++ /dev/null @@ -1,577 +0,0 @@ -;; 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)))))