Update of /project/cl-gsl/cvsroot/cl-gsl In directory common-lisp.net:/tmp/cvs-serv14937 Modified Files: matrix.lisp Log Message: Most functions are now wrapped. Date: Fri Apr 22 04:37:26 2005 Author: edenny Index: cl-gsl/matrix.lisp diff -u cl-gsl/matrix.lisp:1.2 cl-gsl/matrix.lisp:1.3 --- cl-gsl/matrix.lisp:1.2 Mon Apr 18 02:52:16 2005 +++ cl-gsl/matrix.lisp Fri Apr 22 04:37:26 2005 @@ -17,13 +17,30 @@ ;;;; along with this program; if not, write to the Free Software ;;;; Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -(in-package #:cl-gsl-matrix) +(in-package #:cl-gsl-array) + +(defclass gsl-matrix () + ((ptr :accessor ptr :initarg :ptr) + (size-rows :accessor size-rows :initarg :size-rows) + (size-cols :accessor size-cols :initarg :size-cols) + (element-type :accessor element-type :initarg :element-type))) + + +(defclass gsl-matrix-double-float (gsl-matrix) ()) +(defclass gsl-matrix-single-float (gsl-matrix) ()) +(defclass gsl-matrix-integer (gsl-matrix) ()) +(defclass gsl-matrix-complex-double-float (gsl-matrix) ()) +(defclass gsl-matrix-complex-single-float (gsl-matrix) ()) (defmacro def-matrix-type-funcs% (typ) (let ((type-ptr) (type-val) (type-val-ptr) - (type-string)) + (type-string) + (is-real (or (eq typ 'double-float) + (eq typ 'single-float) + (eq typ 'integer)))) + (cond ((eq typ 'double-float) (setq type-ptr 'gsl-matrix-ptr) @@ -55,82 +72,462 @@ `(progn (defun-foreign ,(concatenate 'string "gsl_" type-string "_alloc") - ((size-1 size-t) - (size-2 size-t)) + ((size-rows size-t) + (size-cols size-t)) ,type-ptr) (defun-foreign ,(concatenate 'string "gsl_" type-string "_free") ((m ,type-ptr)) :void) - ))) - -(def-matrix-type-funcs% double-float) -(def-matrix-type-funcs% single-float) -(def-matrix-type-funcs% integer) -(def-matrix-type-funcs% (complex (double-float))) -(def-matrix-type-funcs% (complex (single-float))) - - -;; Function: gsl_matrix * gsl_matrix_alloc (size_t n1, size_t n2) - -;; Function: gsl_matrix * gsl_matrix_calloc (size_t n1, size_t n2) -;; Function: void gsl_matrix_free (gsl_matrix * m) - -;; Function: double gsl_matrix_get (const gsl_matrix * m, size_t i, size_t j) - -;; Function: void gsl_matrix_set (gsl_matrix * m, size_t i, size_t j, double x) - -;; Function: double * gsl_matrix_ptr (gsl_matrix * m, size_t i, size_t j) -;; Function: const double * gsl_matrix_const_ptr (const gsl_matrix * m, size_t i, size_t j) - -;; Function: void gsl_matrix_set_all (gsl_matrix * m, double x) - -;; Function: void gsl_matrix_set_zero (gsl_matrix * m) - -;; Function: void gsl_matrix_set_identity (gsl_matrix * m) - -;; Function: int gsl_matrix_fwrite (FILE * stream, const gsl_matrix * m) - -;; Function: int gsl_matrix_fread (FILE * stream, gsl_matrix * m) - -;; Function: int gsl_matrix_fprintf (FILE * stream, const gsl_matrix * m, const char * format) - -;; Function: int gsl_matrix_fscanf (FILE * stream, gsl_matrix * m) + (defun-foreign ,(concatenate 'string "gsl_" type-string "_memcpy") + ((v1 ,type-ptr) + (v2 ,type-ptr)) + :int) + + (defun-foreign ,(concatenate 'string "gsl_" type-string "_get") + ((m ,type-ptr) + (i size-t) + (j size-t)) + ,type-val) + + (defun-foreign ,(concatenate 'string "gsl_" type-string "_set") + ((m ,type-ptr) + (i size-t) + (j size-t) + (x ,type-val)) + :void) -;; Function: gsl_matrix_view gsl_matrix_submatrix (gsl_matrix * m, size_t k1, size_t k2, size_t n1, size_t n2) -;; Function: gsl_matrix_const_view gsl_matrix_const_submatrix (const gsl_matrix * m, size_t k1, size_t k2, size_t n1, size_t n2) + (defun-foreign ,(concatenate 'string "gsl_" type-string "_set_all") + ((m ,type-ptr) + (x ,type-val)) + :void) -;; Function: gsl_matrix_view gsl_matrix_view_array (double * base, size_t n1, size_t n2) -;; Function: gsl_matrix_const_view gsl_matrix_const_view_array (const double * base, size_t n1, size_t n2) + (defun-foreign ,(concatenate 'string "gsl_" type-string "_set_zero") + ((m ,type-ptr)) + :void) -;; Function: gsl_matrix_view gsl_matrix_view_array_with_tda (double * base, size_t n1, size_t n2, size_t tda) -;; Function: gsl_matrix_const_view gsl_matrix_const_view_array_with_tda (const double * base, size_t n1, size_t n2, size_t tda) + (defun-foreign ,(concatenate 'string "gsl_" type-string "_set_identity") + ((m ,type-ptr)) + :void) -;; Function: gsl_matrix_view gsl_matrix_view_vector (gsl_vector * v, size_t n1, size_t n2) -;; Function: gsl_matrix_const_view gsl_matrix_const_view_vector (const gsl_vector * v, size_t n1, size_t n2) + (defun-foreign ,(concatenate 'string "gsl_" type-string "_isnull") + ((m ,type-ptr)) + :int) -;; Function: gsl_matrix_view gsl_matrix_view_vector_with_tda (gsl_vector * v, size_t n1, size_t n2, size_t tda) -;; Function: gsl_matrix_const_view gsl_matrix_const_view_vector_with_tda (const gsl_vector * v, size_t n1, size_t n2, size_t tda) + (defun-foreign ,(concatenate 'string "wrap_gsl_" type-string "_fwrite") + ((fn :cstring) + (m ,type-ptr)) + :int) + + (defun-foreign ,(concatenate 'string "wrap_gsl_" type-string "_fread") + ((fn :cstring) + (m ,type-ptr)) + :int) + + (defun-foreign ,(concatenate 'string "wrap_gsl_" type-string "_fprintf") + ((fn :cstring) + (m ,type-ptr)) + :int) + + (defun-foreign ,(concatenate 'string "wrap_gsl_" type-string "_fscanf") + ((fn :cstring) + (m ,type-ptr)) + :int) + + (defun-foreign ,(concatenate 'string "gsl_" type-string "_swap") + ((m1 ,type-ptr) + (m2 ,type-ptr)) + :int) + + ,(when is-real + `(progn + (defun-foreign ,(concatenate 'string "gsl_" type-string "_add") + ((ma ,type-ptr) + (mb ,type-ptr)) + :int) + + (defun-foreign ,(concatenate 'string "gsl_" type-string "_sub") + ((ma ,type-ptr) + (mb ,type-ptr)) + :int) + + (defun-foreign ,(concatenate 'string + "gsl_" type-string "_mul_elements") + ((ma ,type-ptr) + (mb ,type-ptr)) + :int) + + (defun-foreign ,(concatenate 'string + "gsl_" type-string "_div_elements") + ((ma ,type-ptr) + (mb ,type-ptr)) + :int) + + (defun-foreign ,(concatenate 'string "gsl_" type-string "_scale") + ((m ,type-ptr) + ;; seems odd that this is :double for all types + (x :double)) + :int) + + (defun-foreign ,(concatenate 'string + "gsl_" type-string "_add_constant") + ((m ,type-ptr) + ;; and again, :double for all types + (x :double)) + :int) + + (defun-foreign ,(concatenate 'string "gsl_" type-string "_max") + ((vec ,type-ptr)) + ,type-val) + + (defun-foreign ,(concatenate 'string "gsl_" type-string "_min") + ((vec ,type-ptr)) + ,type-val) + + + (defun-foreign ,(concatenate 'string + "gsl_" type-string "_max_index") + ((vec ,type-ptr) + (i-ptr size-t-ptr) + (j-ptr size-t-ptr)) + :void) + + (defun-foreign ,(concatenate 'string + "gsl_" type-string "_min_index") + ((vec ,type-ptr) + (i-ptr size-t-ptr) + (j-ptr size-t-ptr)) + :void) + )) + + ,(when (not is-real) + `(progn + (defun-foreign ,(concatenate 'string "gsl_" type-string "_ptr") + ((m ,type-ptr) + (i size-t) + (j size-t)) + (* ,type-val)) + + (defun-foreign ,(concatenate 'string "wrap_gsl_" type-string "_set") + ((m ,type-ptr) + (i size-t) + (j size-t) + (z (* ,type-val))) + :void) + + (defun-foreign ,(concatenate 'string + "wrap_gsl_" type-string "_set_all") + ((m ,type-ptr) + (z (* ,type-val))) + :void) + )) + ))) -;; Function: gsl_vector_view gsl_matrix_row (gsl_matrix * m, size_t i) -;; Function: gsl_vector_const_view gsl_matrix_const_row (const gsl_matrix * m, size_t i) +(def-matrix-type-funcs% double-float) +(def-matrix-type-funcs% single-float) +(def-matrix-type-funcs% integer) +(def-matrix-type-funcs% (complex (double-float))) +(def-matrix-type-funcs% (complex (single-float))) -;; Function: gsl_vector_view gsl_matrix_column (gsl_matrix * m, size_t j) -;; Function: gsl_vector_const_view gsl_matrix_const_column (const gsl_matrix * m, size_t j) -;; Function: gsl_vector_view gsl_matrix_diagonal (gsl_matrix * m) -;; Function: gsl_vector_const_view gsl_matrix_const_diagonal (const gsl_matrix * m) +(defmacro def-matrix-methods% (class-string func-string) + (let ((class-object (kmrcl:concat-symbol "gsl-matrix-" class-string)) + (is-real (or (string= class-string "integer") + (string= class-string "single-float") + (string= class-string "double-float")))) + `(progn -;; Function: gsl_vector_view gsl_matrix_subdiagonal (gsl_matrix * m, size_t k) -;; Function: gsl_vector_const_view gsl_matrix_const_subdiagonal (const gsl_matrix * m, size_t k) + (defmethod alloc ((o ,class-object)) + (setf (ptr o) (,(kmrcl:concat-symbol "gsl-matrix-" func-string "alloc") + (size-rows o) (size-cols o))) + o) + + (defmethod free ((o ,class-object)) + (,(kmrcl:concat-symbol "gsl-matrix-" func-string "free") (ptr o)) + (setf (ptr o) nil) + (setf (size-rows o) nil) + (setf (size-cols o) nil) + (setf (element-type o) nil)) + + (defmethod get-element ((o ,class-object) i &optional j) + (assert (and (typep i 'integer) (>= i 0) (< i (size-rows o)))) + (assert (and (typep j 'integer) (>= j 0) (< j (size-cols o)))) + ,(if is-real + `(,(kmrcl:concat-symbol "gsl-matrix-" func-string "get") + (ptr o) i j) + `(,(kmrcl:concat-symbol "gsl-" func-string ">complex") + (,(kmrcl:concat-symbol "gsl-matrix-" func-string "ptr") + (ptr o) i j)))) + + (defmethod set-element ((o ,class-object) i &optional j x) + (assert (typep x (element-type o))) + (assert (and (typep i 'integer) (>= i 0) (< i (size-rows o)))) + (assert (and (typep j 'integer) (>= j 0) (< j (size-cols o)))) + ,(if is-real + `(,(kmrcl:concat-symbol "gsl-matrix-" func-string "set") + (ptr o) i j x) + `(,(kmrcl:concat-symbol "with-" class-string "->gsl-" func-string + "ptr") (c-ptr x) + (,(kmrcl:concat-symbol "wrap-gsl-matrix-" func-string "set") + (ptr o) i j c-ptr))) + x) + + (defmethod set-all ((o ,class-object) x) + (assert (typep x (element-type o))) + ,(if is-real + `(,(kmrcl:concat-symbol "gsl-matrix-" func-string "set-all") + (ptr o) x) + `(,(kmrcl:concat-symbol "with-" class-string "->gsl-" func-string + "ptr") (c-ptr x) + (,(kmrcl:concat-symbol "wrap-gsl-matrix-" func-string "set-all") + (ptr o) c-ptr))) + o) + + (defmethod set-zero ((o ,class-object)) + (,(kmrcl:concat-symbol "gsl-matrix-" func-string "set-zero") (ptr o)) + o) + + (defmethod set-identity ((o ,class-object)) + (,(kmrcl:concat-symbol "gsl-matrix-" func-string "set-identity") + (ptr o)) + o) + + + (defmethod read-from-binary-file ((o ,class-object) file-name) + (let ((status)) + (uffi:with-cstring (c-file-name file-name) + (setq status + (,(kmrcl:concat-symbol "wrap-gsl-matrix-" func-string + "fread") c-file-name (ptr o)))) + (values o status))) + + (defmethod read-from-file ((o ,class-object) file-name) + (let ((status)) + (uffi:with-cstring (c-file-name file-name) + (setq status + (,(kmrcl:concat-symbol "wrap-gsl-matrix-" func-string + "fscanf") c-file-name (ptr o)))) + (values o status))) + + (defmethod write-to-binary-file (file-name (o ,class-object)) + (let ((status)) + ;; TODO: check if uffi:with-string returns a result, docs unclear. + (uffi:with-cstring (c-file-name file-name) + (setq status + (,(kmrcl:concat-symbol "wrap-gsl-matrix-" func-string + "fwrite") c-file-name (ptr o)))) + status)) + + (defmethod write-to-file (file-name (o ,class-object)) + (let ((status)) + (uffi:with-cstring (c-file-name file-name) + (setq status + (,(kmrcl:concat-symbol "wrap-gsl-matrix-" func-string + "fprintf") c-file-name (ptr o)))) + status)) + + (defmethod swap ((o1 ,class-object) (o2 ,class-object)) + (assert (and (= (size-rows o1) (size-rows o2)) + (= (size-cols o1) (size-cols o2)))) + (let ((status (,(kmrcl:concat-symbol "gsl-matrix-" func-string + "swap") (ptr o1) (ptr o2)))) + (values o1 status))) + + + (defmethod isnull ((o ,class-object)) + (1/0->t/nil (,(kmrcl:concat-symbol "gsl-matrix-" func-string + "isnull") (ptr o)))) + + ,(when is-real + `(progn + (defmethod add ((o1 ,class-object) (o2 ,class-object)) + (assert (and (= (size-rows o1) (size-rows o2)) + (= (size-cols o1) (size-cols o2)))) + (let ((status (,(kmrcl:concat-symbol "gsl-matrix-" func-string + "add") (ptr o1) (ptr o2)))) + (values o1 status))) + + (defmethod sub ((o1 ,class-object) (o2 ,class-object)) + (assert (and (= (size-rows o1) (size-rows o2)) + (= (size-cols o1) (size-cols o2)))) + (let ((status (,(kmrcl:concat-symbol "gsl-matrix-" func-string + "sub") (ptr o1) (ptr o2)))) + (values o1 status))) + + (defmethod mul-elements ((o1 ,class-object) (o2 ,class-object)) + (assert (and (= (size-rows o1) (size-rows o2)) + (= (size-cols o1) (size-cols o2)))) + (let ((status (,(kmrcl:concat-symbol "gsl-matrix-" func-string + "mul-elements") + (ptr o1) (ptr o2)))) + (values o1 status))) + + (defmethod div-elements ((o1 ,class-object) (o2 ,class-object)) + (assert (and (= (size-rows o1) (size-rows o2)) + (= (size-cols o1) (size-cols o2)))) + (let ((status (,(kmrcl:concat-symbol "gsl-matrix-" func-string + "div-elements") + (ptr o1) (ptr o2)))) + (values o1 status))) + + (defmethod scale ((o ,class-object) x) + (assert (typep x (element-type o))) + ;; coerce to double-float looks wrong, but isn't. + (,(kmrcl:concat-symbol "gsl-matrix-" func-string "scale") + (ptr o) (coerce x 'double-float))) + + (defmethod add-constant ((o ,class-object) x) + (assert (typep x (element-type o))) + ;; coerce to double-float looks wrong, but isn't. + (,(kmrcl:concat-symbol "gsl-matrix-" func-string "add-constant") + (ptr o) (coerce x 'double-float))) + + (defmethod max-value ((o ,class-object)) + (,(kmrcl:concat-symbol "gsl-matrix-" func-string "max") (ptr o))) + + (defmethod min-value ((o ,class-object)) + (,(kmrcl:concat-symbol "gsl-matrix-" func-string "min") (ptr o))) + + + (defmethod max-index ((o ,class-object)) + (let ((i-ptr (uffi:allocate-foreign-object 'size-t)) + (j-ptr (uffi:allocate-foreign-object 'size-t))) + (,(kmrcl:concat-symbol "gsl-matrix-" func-string "max-index") + (ptr o) i-ptr j-ptr) + (prog1 + (list (uffi:deref-pointer i-ptr 'size-t) + (uffi:deref-pointer j-ptr 'size-t)) + (uffi:free-foreign-object i-ptr) + (uffi:free-foreign-object j-ptr)))) + + (defmethod min-index ((o ,class-object)) + (let ((i-ptr (uffi:allocate-foreign-object 'size-t)) + (j-ptr (uffi:allocate-foreign-object 'size-t))) + (,(kmrcl:concat-symbol "gsl-matrix-" func-string "min-index") + (ptr o) i-ptr j-ptr) + (prog1 + (list (uffi:deref-pointer i-ptr 'size-t) + (uffi:deref-pointer j-ptr 'size-t)) + (uffi:free-foreign-object i-ptr) + (uffi:free-foreign-object j-ptr)))) + + (defmethod min-max-indicies ((o ,class-object)) + (list (min-index o) (max-index o))) + + (defmethod min-max-values ((o ,class-object)) + (destructuring-bind ((i-min j-min) (i-max j-max)) + (min-max-indicies o) + (list (get-element o i-min j-min) + (get-element o i-max j-max)))) + )) + ))) -;; Function: gsl_vector_view gsl_matrix_superdiagonal (gsl_matrix * m, size_t k) -;; Function: gsl_vector_const_view gsl_matrix_const_superdiagonal (const gsl_matrix * m, size_t k) -;; Function: int gsl_matrix_memcpy (gsl_matrix * dest, const gsl_matrix * src) +(def-matrix-methods% "integer" "int-") +(def-matrix-methods% "single-float" "float-") +(def-matrix-methods% "double-float" "") +(def-matrix-methods% "complex-single-float" "complex-float-") +(def-matrix-methods% "complex-double-float" "complex-") + + +(defun make-matrix (size-rows size-cols + &key (element-type 'double-float) initial-element + initial-contents from-file from-binary-file) + (assert (and (typep size-rows 'integer) (> size-rows 0) )) + (assert (and (typep size-cols 'integer) (> size-cols 0) )) + (assert (find element-type '(integer single-float double-float + (complex (single-float)) + (complex (double-float))) :test #'equal)) + (let ((m (cond + ((eq element-type 'integer) + (make-instance 'gsl-matrix-integer + :size-rows size-rows :size-cols size-cols + :element-type element-type)) + ((eq element-type 'double-float) + (make-instance 'gsl-matrix-double-float + :size-rows size-rows :size-cols size-cols + :element-type element-type)) + ((eq element-type 'single-float) + (make-instance 'gsl-matrix-single-float + :size-rows size-rows :size-cols size-cols + :element-type element-type)) + ((equal element-type '(complex (double-float))) + (make-instance 'gsl-matrix-complex-double-float + :size-rows size-rows :size-cols size-cols + :element-type element-type)) + ((equal element-type '(complex (single-float))) + (make-instance 'gsl-matrix-complex-single-float + :size-rows size-rows :size-cols size-cols + :element-type element-type)) + (t + (error "should never get here."))))) + (alloc m) + (cond + ((and initial-element initial-contents from-file from-binary-file) + (error "can only define one of the keys: initial-element, initial-contents, from-file, from-binary-file.")) + (initial-element + (set-all m initial-element)) + (initial-contents + (cond + ((arrayp initial-contents) + (dotimes (i size-rows) + (dotimes (j size-cols) + (set-element m i j (aref initial-contents i j))))) + (t + (error "initial-contents must be an array.")))) + (from-file + (read-from-file m from-file)) + (from-binary-file + (read-from-binary-file m from-binary-file))) + m)) + + +(defmacro with-matrix + ((m size-rows size-cols &key element-type initial-element initial-contents + from-file from-binary-file) &body body) + `(let ((,m (make-matrix ,size-rows ,size-cols + :element-type (or ,element-type 'double-float) + :initial-element ,initial-element + :initial-contents ,initial-contents + :from-file ,from-file + :from-binary-file ,from-binary-file))) + (unwind-protect + (progn ,@body) + (free ,m)))) + + +(defmacro def-matrix-copy-method% (class-string func-string) + (let ((class-object (kmrcl:concat-symbol "gsl-matrix-" class-string))) + `(defmethod copy ((o ,class-object)) + (let* ((o-copy (make-matrix (size-rows o) (size-cols o) + :element-type (element-type o))) + (status (,(kmrcl:concat-symbol "gsl-matrix-" func-string + "memcpy") (ptr o-copy) (ptr o)))) + (values o-copy status))))) + +(def-matrix-copy-method% "integer" "int-") +(def-matrix-copy-method% "single-float" "float-") +(def-matrix-copy-method% "double-float" "") +(def-matrix-copy-method% "complex-single-float" "complex-float-") +(def-matrix-copy-method% "complex-double-float" "complex-") + + +(defmacro with-matrix-copy ((m-dest m-src) &body body) + `(let ((,m-dest (copy ,m-src))) + (unwind-protect + ,@body + (free ,m-dest)))) + + +(defun gsl-matrix->lisp-array (m) + (let ((a (make-array (list (size-rows m) (size-cols m)) + :element-type (element-type m)))) + (dotimes (i (size-rows m)) + (dotimes (j (size-cols m)) + (setf (aref a i j) (get-element m i j)))) + a)) -;; Function: int gsl_matrix_swap (gsl_matrix * m1, gsl_matrix * m2) ;; Function: int gsl_matrix_get_row (gsl_vector * v, const gsl_matrix * m, size_t i) @@ -150,30 +547,3 @@ ;; Function: int gsl_matrix_transpose_memcpy (gsl_matrix * dest, const gsl_matrix * src) ;; Function: int gsl_matrix_transpose (gsl_matrix * m) - -;; Function: int gsl_matrix_add (gsl_matrix * a, const gsl_matrix * b) - -;; Function: int gsl_matrix_sub (gsl_matrix * a, const gsl_matrix * b) - -;; Function: int gsl_matrix_mul_elements (gsl_matrix * a, const gsl_matrix * b) - -;; Function: int gsl_matrix_div_elements (gsl_matrix * a, const gsl_matrix * b) - -;; Function: int gsl_matrix_scale (gsl_matrix * a, const double x) - -;; Function: int gsl_matrix_add_constant (gsl_matrix * a, const double x) - -;; Function: double gsl_matrix_max (const gsl_matrix * m) - -;; Function: double gsl_matrix_min (const gsl_matrix * m) - -;; Function: void gsl_matrix_minmax (const gsl_matrix * m, double * min_out, double * max_out) - -;; Function: void gsl_matrix_max_index (const gsl_matrix * m, size_t * imax, size_t * jmax) - -;; Function: void gsl_matrix_min_index (const gsl_matrix * m, size_t * imax, size_t * jmax) - -;; Function: void gsl_matrix_minmax_index (const gsl_matrix * m, size_t * imin, size_t * imax) - -;; Function: int gsl_matrix_isnull (const gsl_matrix * m) -