Menu
▾
▴
[Matlisp-commit] [matlisp-git]matlisp branch matlisp-cffi updated. 2012-02-23-71-gff26318
|
From: Akshay S. <ak...@us...> - 2012-03-24 08:32:04
|
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "matlisp".
The branch, matlisp-cffi has been updated
via ff263186ffc1a8443f5733cc975ba2e7c66d2206 (commit)
via f53e544eff8af4aa8fdd302a1adc98fed1b5aa35 (commit)
from 9bfeec0a8b2e5604b2ce6b7ad6be62c3fd3f09c4 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit ff263186ffc1a8443f5733cc975ba2e7c66d2206
Author: Akshay Srinivasan <aks...@gm...>
Date: Sat Mar 24 13:57:22 2012 +0530
o Implemented support for callbacks.
o Stated to using new protocol to append "~" to functions
which return matrices which share the store.
o Lots of tweaks to the FFI.
diff --git a/matlisp.asd b/matlisp.asd
index 2529cd7..2c8019c 100644
--- a/matlisp.asd
+++ b/matlisp.asd
@@ -147,8 +147,9 @@
(:file "help")
(:file "special")
(:file "reader")
- ;;(:file "trans")
- ;;(:file "realimag")
+ (:file "trans")
+ (:file "realimag")
+ (:file "submat")
(:file "reshape")
(:file "join")
(:file "svd")
diff --git a/packages.lisp b/packages.lisp
index e4fe2da..94f67e0 100644
--- a/packages.lisp
+++ b/packages.lisp
@@ -160,22 +160,27 @@
#:zip-eq
#:cut-cons-chain!
#:when-let
+ #:if-let
#:if-ret
#:get-arg
#:nconsc
#:with-gensyms
#:slot-values
- #:mlet*))
+ #:mlet*
+ #:recursive-append
+ ;;
+ #:foreign-vector #:make-foreign-vector #:foreign-vector-p
+ #:fv-ref #:fv-pointer #:fv-size #:fv-type))
(defpackage :fortran-ffi-accessors
+ (:nicknames :ffi)
#+:cmu (:use :common-lisp :c-call :cffi :utilities)
#+:sbcl (:use :common-lisp :sb-alien :sb-c :cffi :utilities)
#+:allegro (:use :common-lisp :cffi :utilities)
#+(not (or sbcl cmu allegro)) (:use :common-lisp :cffi :utilities)
(:export
;; interface functions
- #:def-fortran-routine
- #:incf-sap
+ #:def-fortran-routine
#:with-vector-data-addresses
)
(:documentation "Fortran foreign function interface"))
@@ -315,14 +320,22 @@
#:store #:store-size
;;Generic functions on standard-matrix
#:fill-matrix
- #:ctranspose! #:ctranspose #:transpose #:transpose!
#:row-or-col-vector-p #:row-vector-p #:col-vector-p
- #:row #:col #:diag #:sub-matrix
+ ;;Submatrix ops
+ #:row~ #:row
+ #:col~ #:col
+ #:diag~ #:diag
+ #:sub-matrix~ #:sub-matrix
+ ;;Transpose
+ #:transpose~ #:transpose! #:transpose
+ #:ctranspose! #:ctranspose
;;Real-double-matrix
#:real-matrix #:real-matrix-element-type #:real-matrix-store-type
;;Complex-double-matrix
#:complex-matrix #:complex-matrix-element-type #:complex-matrix-store-type #:complex-coerce #:complex-double-float
- #:mrealpart #:mimagpart
+ ;;Real and imaginary parts
+ #:mrealpart~ #:mrealpart #:real
+ #:mimagpart~ #:mimagpart #:imag
;;
"CONVERT-TO-LISP-ARRAY"
"DOT"
@@ -399,7 +412,6 @@
"POTRF!"
"POTRS!"
"RAND"
- "REAL"
"RESHAPE!"
"RESHAPE"
"SAVE-MATLISP"
diff --git a/src/axpy.lisp b/src/axpy.lisp
index 662cfd6..0891c60 100644
--- a/src/axpy.lisp
+++ b/src/axpy.lisp
@@ -134,11 +134,11 @@ don't know how to coerce COMPLEX to REAL"))
(generate-typed-axpy!-func complex-double-axpy!-typed complex-double-float complex-matrix-store-type complex-matrix blas:zaxpy)
(defmethod axpy! ((alpha cl:real) (x real-matrix) (y complex-matrix))
- (real-double-axpy!-typed (coerce alpha 'double-float) x (mrealpart y)))
+ (real-double-axpy!-typed (coerce alpha 'double-float) x (mrealpart~ y)))
(defmethod axpy! ((alpha complex) (x real-matrix) (y complex-matrix))
- (real-double-axpy!-typed (coerce (realpart alpha) 'double-float) x (mrealpart y))
- (real-double-axpy!-typed (coerce (imagpart alpha) 'double-float) x (mimagpart y)))
+ (real-double-axpy!-typed (coerce (realpart alpha) 'double-float) x (mrealpart~ y))
+ (real-double-axpy!-typed (coerce (imagpart alpha) 'double-float) x (mimagpart~ y)))
(defmethod axpy! ((alpha number) (x complex-matrix) (y complex-matrix))
(complex-double-axpy!-typed (complex-coerce alpha) x y))
diff --git a/src/complex-matrix.lisp b/src/complex-matrix.lisp
index 3ad8280..928d43f 100644
--- a/src/complex-matrix.lisp
+++ b/src/complex-matrix.lisp
@@ -69,75 +69,6 @@
(aref store (+ 1 (* 2 idx))) (imagpart coerced-value))))
;;
-(defmethod transpose ((matrix complex-matrix))
- (mlet* (((hd nr nc rs cs st) (slot-values matrix '(head number-of-rows number-of-cols row-stride col-stride store))
- :type (fixnum fixnum fixnum fixnum fixnum (complex-matrix-store-type *))))
- (make-instance 'sub-complex-matrix
- :nrows nc :ncols nr
- :store st
- :head hd
- :row-stride cs :col-stride rs
- :parent matrix)))
-
-;;
-(defmethod sub-matrix ((matrix complex-matrix) (origin list) (dim list))
- (destructuring-bind (o-i o-j) origin
- (destructuring-bind (nr-s nc-s) dim
- (mlet* (((hd nr nc rs cs st) (slot-values matrix '(head number-of-rows number-of-cols row-stride col-stride store))
- :type (fixnum fixnum fixnum fixnum fixnum (complex-matrix-store-type *))))
- (unless (and (< -1 o-i (+ o-j nr-s) nr) (< -1 o-j (+ o-j nc-s) nc))
- (error "Bad index and/or size.
-Cannot create a sub-matrix of size (~a ~a) starting at (~a ~a)" nr-s nc-s o-i o-j))
- (make-instance 'sub-complex-matrix
- :nrows nr-s :ncols nc-s
- :store st
- :head (store-indexing o-i o-j hd rs cs)
- :row-stride rs :col-stride cs)))))
-
-;;
-(defmethod row ((matrix complex-matrix) (i fixnum))
- (mlet* (((hd nr nc rs cs st) (slot-values matrix '(head number-of-rows number-of-cols row-stride col-stride store))
- :type (fixnum fixnum fixnum fixnum fixnum (complex-matrix-store-type *))))
- (unless (< -1 i nr)
- (error "Index ~a is outside the valid range for the given matrix." i))
- (make-instance 'sub-complex-matrix
- :nrows 1 :ncols nc
- :store st
- :head (store-indexing i 0 hd rs cs)
- :row-stride rs :col-stride cs)))
-
-;;
-(defmethod col ((matrix complex-matrix) (j fixnum))
- (mlet* (((hd nr nc rs cs st) (slot-values matrix '(head number-of-rows number-of-cols row-stride col-stride store))
- :type (fixnum fixnum fixnum fixnum fixnum (complex-matrix-store-type *))))
- (unless (< -1 j nc)
- (error "Index ~a is outside the valid range for the given matrix." j))
- (make-instance 'sub-complex-matrix
- :nrows nr :ncols 1
- :store st
- :head (store-indexing 0 j hd rs cs)
- :row-stride rs :col-stride cs)))
-
-;;
-(defmethod diag ((matrix complex-matrix) &optional (d 0))
- (declare (type fixnum d))
- (mlet* (((hd nr nc rs cs st) (slot-values matrix '(head number-of-rows number-of-cols row-stride col-stride store))
- :type (fixnum fixnum fixnum fixnum fixnum (complex-matrix-store-type *)))
- ((f-i f-j) (if (< d 0)
- (values (- d) 0)
- (values 0 d))
- :type (fixnum fixnum)))
- (unless (and (< -1 f-i nr) (< -1 f-j nc))
- (error "Index ~a is outside the valid range for the given matrix." d))
- (let ((d-s (min (- nr f-i) (- nc f-j))))
- (declare (type fixnum d-s))
- (make-instance 'sub-complex-matrix
- :nrows 1 :ncols d-s
- :store st
- :head (store-indexing f-i f-j hd rs cs)
- :row-stride 1 :col-stride (+ rs cs)))))
-
-;;
(declaim (inline allocate-complex-store))
(defun allocate-complex-store (size)
(make-array (* 2 size) :element-type 'complex-matrix-element-type
@@ -357,31 +288,3 @@ Cannot create a sub-matrix of size (~a ~a) starting at (~a ~a)" nr-s nc-s o-i o-
(make-complex-matrix-dim n m)))
(t
(error "require 1 or 2 arguments to make a matrix")))))
-
-;;
-
-(defun mrealpart (mat)
- (typecase mat
- (real-matrix mat)
- (complex-matrix (make-instance 'sub-real-matrix
- :parent mat :store (store mat)
- :nrows (nrows mat) :ncols (ncols mat)
- :row-stride (* 2 (row-stride mat)) :col-stride (* 2 (col-stride mat))
- :head (* 2 (head mat))))
- (number (cl:realpart mat))))
-
-(defun mimagpart (mat)
- (typecase mat
- (real-matrix nil)
- (complex-matrix (make-instance 'sub-real-matrix
- :parent mat :store (store mat)
- :nrows (nrows mat) :ncols (ncols mat)
- :row-stride (* 2 (row-stride mat)) :col-stride (* 2 (col-stride mat))
- :head (+ 1 (* 2 (head mat)))))
- (number (cl:imagpart mat))))
-
-(defun mconjugate! (mat)
- (typecase mat
- (real-matrix mat)
- (complex-matrix (scal! -1d0 (mimagpart mat))))
- mat)
\ No newline at end of file
diff --git a/src/copy.lisp b/src/copy.lisp
index 344999a..aa3cbda 100644
--- a/src/copy.lisp
+++ b/src/copy.lisp
@@ -219,8 +219,8 @@ don't know how to coerce a COMPLEX to a REAL"))
(complex-double-copy!-typed x y))
(defmethod copy! ((x real-matrix) (y complex-matrix))
- (real-double-copy!-typed x (mrealpart y))
- (scal! 0d0 (mimagpart y))
+ (real-double-copy!-typed x (mrealpart~ y))
+ (scal! 0d0 (mimagpart~ y))
y)
(defmethod copy! ((x number) (y complex-matrix))
diff --git a/src/dlsode.lisp b/src/dlsode.lisp
index 3883519..a284d50 100644
--- a/src/dlsode.lisp
+++ b/src/dlsode.lisp
@@ -1,15 +1,7 @@
-(in-package "MATLISP")
-#+nil
-(progn
-(asdf:oos 'asdf:load-op :cffi)
-
-(load "f77-mangling.lisp")
-(load "cffi-helpers.lisp")
-(load "ffi-cffi.lisp")
-)
+(in-package #:matlisp)
(cffi:define-foreign-library libodepack
- (:unix #.(translate-logical-pathname
+ #+nil(:unix #.(translate-logical-pathname
(merge-pathnames "matlisp:lib;libodepack"
*shared-library-pathname-extension*)))
(t (:default "libodepack")))
@@ -17,13 +9,23 @@
(cffi:use-foreign-library libodepack)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+#+nil(def-fortran-routine testde :void
+ (field (:callback :void
+ (c-neq :integer :input)
+ (c-t :double-float :input)
+ (c-y (* :double-float :size c-neq) :input)
+ (c-ydot (* :double-float :size c-neq) :output)))
+ (neq :integer :input)
+ (y (* :double-float) :input-output))
+
+
(def-fortran-routine dlsode :void
"DLSODE in ODEPACK"
(field (:callback :void
(c-neq :integer :input)
(c-t :double-float :input)
- (c-y (* :double-float) :input)
- (c-ydot (* :double-float) :output)))
+ (c-y (* :double-float :size c-neq) :input)
+ (c-ydot (* :double-float :size c-neq) :output)))
(neq :integer :input)
(y (* :double-float) :input-output)
(ts :double-float :input-output)
@@ -41,31 +43,15 @@
(jacobian (:callback :void
(c-neq :integer :input)
(c-t :double-float :input)
- (c-y (* :double-float) :input)
+ (c-y (* :double-float :size c-neq) :input)
(c-upper-bandwidth :integer :input)
(c-lower-bandwidth :integer :input)
- (c-pd (* :double-float) :output)
+ (c-pd (* :double-float :size (* c-neq c-neq)) :output)
(c-nrowpd :integer :input)))
(mf :integer :input))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun lsode-evolve (field y t-array report)
- ;; Use gensym ? Will have to use a macrolet.
- (cffi:defcallback *evolve-callback* :void ((c-neq :pointer :int)
- (c-tc :pointer :double)
- (c-y :pointer :double)
- (c-ydot :pointer :double))
- (let* ((neq (cffi:mem-aref c-neq :int))
- (y (make-array neq :element-type 'double-float :initial-element 0d0))
- (ts (cffi:mem-aref c-tc :double)))
- ;; Copy things to simple-arrays
- (loop for i from 0 below neq
- do (setf (aref y i) (cffi:mem-aref c-y :double i)))
- ;; Assume form of field
- (let ((ydot (funcall field ts y)))
- ;; Copy ydot back
- (loop for i from 0 below neq
- do (setf (cffi:mem-aref c-ydot :double i) (aref ydot i))))))
;;
(let* ((neq (length y))
(lrw (+ 22 (* 9 neq) (* neq neq) 5))
@@ -86,21 +72,22 @@
do (progn
(setq tout (aref t-array i))
(multiple-value-bind (y-out ts-out istate-out rwork-out iwork-out)
- (dlsode (cffi:callback *evolve-callback*) neq y ts tout itol rtol atol itask istate iopt rwork lrw iwork liw (cffi:null-pointer) mf)
+ (dlsode field neq y ts tout itol rtol atol itask istate iopt rwork lrw iwork liw (cffi:null-pointer) mf)
(setq ts ts-out)
(setq istate istate-out))
(funcall report ts y)))))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-(defun pend-field (ts y)
- (make-array 2 :element-type 'double-float :initial-contents `(,(aref y 1) ,(- (sin (aref y 0))))))
+(defun pend-field (neq time y ydot)
+ (setf (fv-ref ydot 0) (fv-ref y 1)
+ (fv-ref ydot 1) (- (sin (fv-ref y 0)))))
(defun pend-report (ts y)
(format t "~A ~A ~A ~%" ts (aref y 0) (aref y 1)))
(defvar y (make-array 2 :element-type 'double-float :initial-contents `(,(/ pi 2) 0d0)))
-(lsode-evolve #'pend-field y #(0d0 1d0) #'pend-report)
+;; (lsode-evolve #'pend-field y #(0d0 1d0 2d0) #'pend-report)
;; Should return
;; 1.0d0 1.074911802207049d0 -0.975509986605856d0
;; 2.0d0 -0.20563950412081608d0 -1.3992359518735706d0
diff --git a/src/dot.lisp b/src/dot.lisp
index badff89..5a5d896 100644
--- a/src/dot.lisp
+++ b/src/dot.lisp
@@ -139,7 +139,7 @@
(let ((realpart (ddot nxm store-x 1 store-y 2))
(imagpart (with-vector-data-addresses ((addr-x store-x)
(addr-y store-y))
- (incf-sap :double-float addr-y)
+ (incf-sap addr-y :double-float)
(ddot nxm addr-x 1 addr-y 2))))
(declare (type complex-matrix-element-type realpart imagpart))
@@ -192,7 +192,7 @@
(let ((realpart (ddot nxm store-x 2 store-y 1))
(imagpart (with-vector-data-addresses ((addr-x store-x)
(addr-y store-y))
- (incf-sap :double-float addr-x)
+ (incf-sap addr-x :double-float)
(ddot nxm addr-x 2 addr-y 1))))
(declare (type complex-matrix-element-type realpart imagpart))
diff --git a/src/ffi-cffi-interpreter-specific.lisp b/src/ffi-cffi-interpreter-specific.lisp
index 75afa9a..eda8343 100644
--- a/src/ffi-cffi-interpreter-specific.lisp
+++ b/src/ffi-cffi-interpreter-specific.lisp
@@ -82,14 +82,16 @@ Returns
#+(or sbcl cmu ccl)
(defmacro with-vector-data-addresses (vlist &body body)
- "WITH-VECTOR-DATA-ADDRESSES (var-list &body body)
+"
+ WITH-VECTOR-DATA-ADDRESSES (var-list &body body)
- Execute the body with the variables in VAR-LIST appropriately bound.
- VAR-LIST should be a list of pairs. The first element is the address
- of the desired object; the second element is the variable whose address
- we want.
+ Execute the body with the variables in VAR-LIST appropriately bound.
+ VAR-LIST should be a list of pairs. The first element is the address
+ of the desired object; the second element is the variable whose address
+ we want.
- Garbage collection is also disabled while executing the body."
+ Garbage collection is also disabled while executing the body.
+"
;; We wrap everything inside a WITHOUT-GCING form to inhibit garbage
;; collection to avoid complications that may arise during a
;; collection while in a fortran call.
@@ -103,13 +105,7 @@ Returns
(let (,@(mapcar #'(lambda (lst)
(destructuring-bind (addr-var var &key inc-type inc) lst
`(,addr-var ,@(if inc
- `((cffi:inc-pointer (vector-data-address ,var)
- ,@(case inc-type
- (:double-float `((* ,inc 8)))
- (:single-float `((* ,inc 4)))
- (:complex-double-float `((* ,inc 16)))
- (:complex-single-float `((* ,inc 8)))
- (t `(,inc)))))
+ `((inc-sap (vector-data-address ,var) ,inc-type ,inc))
`((vector-data-address ,var))))))
vlist))
,@body))))
diff --git a/src/ffi-cffi.lisp b/src/ffi-cffi.lisp
index b1d8bcc..5146374 100644
--- a/src/ffi-cffi.lisp
+++ b/src/ffi-cffi.lisp
@@ -11,97 +11,102 @@
(in-package "FORTRAN-FFI-ACCESSORS")
-#+(or)
(defconstant +ffi-types+ '(:single-float :double-float
:complex-single-float :complex-double-float
:integer :long
:string
:callback))
-
(defconstant +ffi-styles+ '(:input :input-value :workspace
;;
:input-output :output :workspace-output))
+
;; Create objects on the heap and run some stuff.
-(defmacro with-foreign-objects-heap-ed (declarations &rest body)
- "
-Allocate \"objects\" on the heap and run the \"body\" of code.
+(defmacro with-foreign-objects-heaped (declarations &rest body)
+"
+ Allocate \"objects\" on the heap and run the \"body\" of code.
-with-foreign-objects-heap-ed (declarations) &rest body
-binding := {(var type &optional count &key (initial-contents nil))}*
+ with-foreign-objects-heap-ed (declarations) &rest body
+ binding := {(var type &optional count &key (initial-contents nil))}*
-Example:
->> (with-foreign-objects-heap-ed ((x :int :count 10 :initial-element 2))
- (+ (cffi:mem-aref x :int 2) 1))
-3
+ Example:
+ >> (with-foreign-objects-heaped ((x :int :count 10 :initial-element 2))
+ (+ (cffi:mem-aref x :int 2) 1))
+ 3
>>
- "
- (let ((ret (gensym)))
- ;; Allocate objects from the heap
- `(let* (,@(mapcar (lambda (decl) (list (car decl) `(cffi:foreign-alloc ,@(cdr decl))))
- declarations)
- ;; Store result temporarily
- (,ret (progn ,@body)))
- ;;Free C objects
- ,@(mapcar (lambda (decl) `(cffi:foreign-free ,(car decl)))
- declarations)
- ,ret)))
+"
+;; Allocate objects from the heap
+ (recursive-append
+ (when declarations
+ `(let (,@(mapcar (lambda (decl) (let ((var (car decl)))
+ (check-type var symbol)
+ `(,var (cffi:foreign-alloc ,@(cdr decl)))))
+ declarations))))
+ ;; Store result and then free foreign-objects
+ (when declarations
+ `(multiple-value-prog1))
+ `((progn
+ ,@body)
+ ;;Free C objects
+ ,@(mapcar (lambda (decl) `(cffi:foreign-free ,(car decl)))
+ declarations))))
;; Create objects on the stack and run the "body" of code.
-(defmacro with-foreign-objects-stack-ed (declarations &rest body)
- "
-Allocate \"objects\" on the stack and run the \"body\" of code.
-
-with-foreign-objects-stack-ed (declarations) &rest body
-binding := {(var type &optional count &key (initial-contents nil))}*
-
-Example:
->> (with-foreign-objects-stack-ed ((x :int :count 10 :initial-element 2))
- (+ (cffi:mem-aref x :int 2) 1))
-3
->>
- "
- (if (null declarations)
- `(progn ,@body)
- (let ((wfo-decl nil)
- (wfo-body nil)
- (wfo-before nil))
- (loop for decl in declarations
- do (destructuring-bind (var type &key (count 1) initial-element initial-contents) decl
- ;;Make sure the var and type are symbols;;
- (check-type var symbol)
- (check-type type symbol)
- (when (and initial-element initial-contents)
- (error "Cannot apply both :initial-element and :initial-contents at the same time."))
- ;;
- (if (eq count 1)
- (progn
- ;; Count defaults to one in with-foreign-objects
- (nconsc wfo-decl `((,var ,type)))
- (if (or initial-element initial-contents)
- (nconsc wfo-body `((setf (cffi:mem-aref ,var ,type 0) ,@(cond
- (initial-element `(,initial-element))
- (initial-contents `((car ,initial-contents)))))))))
- ;;
- (let ((decl-count (gensym))
- (decl-init (gensym))
- (loop-var (gensym)))
- ;;
- (nconsc wfo-before `((,decl-count ,count)))
- (nconsc wfo-before `((,decl-init ,(or initial-element initial-contents))))
- ;;
- (nconsc wfo-decl `((,var ,type ,decl-count)))
- (if (or initial-element initial-contents)
- (nconsc wfo-body `((loop for ,loop-var from 0 below ,decl-count
- do (setf (cffi:mem-aref ,var ,type ,loop-var) ,@(cond
- (initial-element `(,decl-init))
-
- (initial-contents `((elt ,decl-init ,loop-var)))))))))))))
- `(let (,@wfo-before)
- (cffi:with-foreign-objects (,@wfo-decl)
- ,@wfo-body
- ,@body)))))
+(defmacro with-foreign-objects-stacked (declarations &rest body)
+"
+ Allocate \"objects\" on the stack and run the \"body\" of code.
+
+ with-foreign-objects-stacked (declarations) &rest body
+ binding := {(var type &optional count &key (initial-contents nil))}*
+
+ Example:
+ >> (with-foreign-objects-stacked ((x :int :count 10 :initial-element 2))
+ (+ (cffi:mem-aref x :int 2) 1))
+ 3
+ >>
+"
+ (let ((wfo-decl nil)
+ (wfo-body nil)
+ (wfo-before nil))
+ (dolist (decl declarations)
+ (destructuring-bind (var type &key (count 1) initial-element initial-contents) decl
+ ;;Make sure the var and type are symbols;;
+ (check-type var symbol)
+ (check-type type symbol)
+ (when (and initial-element initial-contents)
+ (error "Cannot apply both :initial-element and :initial-contents at the same time."))
+ ;;
+ (if (eq count 1)
+ (progn
+ ;; Count defaults to one in with-foreign-objects
+ (nconsc wfo-decl `((,var ,type)))
+ (if (or initial-element initial-contents)
+ (nconsc wfo-body `((setf (cffi:mem-aref ,var ,type 0) ,@(cond
+ (initial-element `(,initial-element))
+ (initial-contents `((elt ,initial-contents 0)))))))))
+ ;;
+ (let ((decl-count (gensym))
+ (decl-init (gensym))
+ (loop-var (gensym)))
+ ;;
+ (nconsc wfo-before `((,decl-count ,count)))
+ (nconsc wfo-before `((,decl-init ,(or initial-element initial-contents))))
+ ;;
+ (nconsc wfo-decl `((,var ,type ,decl-count)))
+ (if (or initial-element initial-contents)
+ (nconsc wfo-body `((dotimes (,loop-var ,decl-count)
+ (setf (cffi:mem-aref ,var ,type ,loop-var) ,@(cond
+ (initial-element `(,decl-init))
+ (initial-contents `((elt ,decl-init ,loop-var)))))))))))))
+ (recursive-append
+ (when wfo-before
+ `(let (,@wfo-before)))
+ (if wfo-decl
+ `(cffi:with-foreign-objects (,@wfo-decl))
+ `(progn))
+ `(,@wfo-body
+ ,@body))))
;; Get the equivalent CFFI type.
;; If the type is an array, get the type of the array element type.
@@ -201,8 +206,7 @@ Example:
(:string
;; String lengths are appended to the function arguments,
;; passed by value.
- (pushnew `(,(scat "LEN-" name) ,@(->cffi-type :integer))
- aux-pars)
+ (nconsc aux-pars `((,(scat "LEN-" name) ,@(->cffi-type :integer))))
`(,name ,@(->cffi-type :string)))
(t
`(,name ,@(get-read-in-type type style))))))
@@ -214,14 +218,14 @@ Example:
;; Call defcfun to define the foreign function.
;; Also creates a nice lisp helper function.
(defmacro def-fortran-routine (func-name return-type &rest body)
- (multiple-value-bind (name fortran-name) (if (listp func-name)
- (values (cadr func-name) (car func-name))
- (values func-name (make-fortran-name func-name)))
+ (multiple-value-bind (fortran-name name) (if (listp func-name)
+ (values (car func-name) (cadr func-name))
+ (values (make-fortran-name func-name) func-name))
(let* ((lisp-name (make-fortran-ffi-name `,name))
(hack-return-type `,return-type)
(hack-body `(,@body))
(hidden-var-name nil))
-
+ ;;
(multiple-value-bind (doc pars)
(parse-doc-&-parameters `(,@body))
(when (member hack-return-type '(:complex-single-float :complex-double-float))
@@ -236,7 +240,7 @@ Example:
(,hidden-var-name ,hack-return-type :output)
,@pars))
(setq hack-return-type :void)))
-
+
`(eval-when (load eval compile)
(progn
;; Removing 'inlines' It seems that CMUCL has a problem with
@@ -257,6 +261,7 @@ Example:
(return-vars nil)
(array-vars nil)
(ref-vars nil)
+ (callback-code nil)
;;
(defun-args nil)
(defun-keyword-args nil)
@@ -265,133 +270,230 @@ Example:
;;
(ffi-args nil)
(aux-ffi-args nil))
- (loop for decl in pars
- do (destructuring-bind (var type &optional style) decl
- (let ((ffi-var nil)
- (aux-var nil))
- (cond
- ;; Callbacks are tricky because the data inside
- ;; pointer arrays will need to be copied without
- ;; implicit knowledge of the size of the array.
- ;; This is usually taken care of by special data
- ;; structure - ala GSL - or by passing additional
- ;; arguments to the callback to apprise it of the
- ;; bounds on the arrays.
- ;; TODO: Add support for declaring array dimensions
- ;; in the callback declaration.
- ((callback-type-p type)
- (setq ffi-var var))
- ;; Can't really enforce "style" when given an array.
- ;; Complex numbers do not latch onto this case, they
- ;; are passed by value.
- ((array-p type)
- (setq ffi-var (scat "ADDR-" var))
- (nconsc array-vars `((,ffi-var ,var)))
- ;;
- (when-let (arg (get-arg :inc type))
- (nconsc defun-keyword-args
- `((,arg 0)))
- (nconc (car (last array-vars)) `(:inc-type ,(cadr type) :inc ,arg))))
- ;; Strings
- ((string-p type)
- (setq ffi-var var)
- (setq aux-var (scat "LEN-" var))
- (nconsc aux-args `((,aux-var (length (the string ,var))))))
- ;; Pass-by-value variables
- ((eq style :input-value)
- (setq ffi-var var))
- ;; Pass-by-reference variables
- (t
- (cond
- ;; Makes more sense to copy complex numbers into
- ;; arrays, rather than twiddling around with lisp
- ;; memory internals.
- ((member type '(:complex-single-float :complex-double-float))
- (setq ffi-var (scat "ADDR-REAL-CAST-" var))
- (nconsc ref-vars
- `((,ffi-var ,(second (->cffi-type type)) :count 2 :initial-contents (list (realpart ,var) (imagpart ,var))))))
- (t
- (setq ffi-var (scat "REF-" var))
- (nconsc ref-vars
- `((,ffi-var ,@(->cffi-type type) :initial-element ,var)))))))
- ;; Output variables
- (when (and (output-p style) (not (eq type :string)))
- (nconsc return-vars
- `((,ffi-var ,var ,type))))
- ;; Arguments for the lisp wrapper
- (when (not (eq var hidden-var-name))
- (nconsc defun-args
- `(,var)))
- ;; Arguments for the FFI function
- (nconsc ffi-args
- `(,ffi-var))
- ;; Auxillary arguments for FFI
- (when (not (null aux-var))
- (nconsc aux-ffi-args
- `(,aux-var))))))
- ;;Return the function definition
+ (dolist (decl pars)
+ (destructuring-bind (var type &optional style) decl
+ (let ((ffi-var nil)
+ (aux-var nil))
+ (cond
+ ;; Callbacks are tricky.
+ ((callback-type-p type)
+ (let* ((callback-name (gensym (symbol-name var)))
+ (c-callback-code (def-fortran-callback var callback-name (second type) (cddr type))))
+ (nconsc callback-code c-callback-code)
+ (setq ffi-var `(cffi:callback ,callback-name))))
+ ;; Can't really enforce "style" when given an array.
+ ;; Complex numbers do not latch onto this case, they
+ ;; are passed by value.
+ ((array-p type)
+ (setq ffi-var (scat "ADDR-" var))
+ (nconsc array-vars `((,ffi-var ,var)))
+ ;;
+ (when-let (arg (get-arg :inc type))
+ (nconsc defun-keyword-args
+ `((,arg 0)))
+ (nconc (car (last array-vars)) `(:inc-type ,(cadr type) :inc ,arg))))
+ ;; Strings
+ ((string-p type)
+ (setq ffi-var var)
+ (setq aux-var (scat "LEN-" var))
+ (nconsc aux-args `((,aux-var (length (the string ,var))))))
+ ;; Pass-by-value variables
+ ((eq style :input-value)
+ (setq ffi-var var))
+ ;; Pass-by-reference variables
+ (t
+ (cond
+ ;; Makes more sense to copy complex numbers into
+ ;; arrays, rather than twiddling around with lisp
+ ;; memory internals.
+ ((member type '(:complex-single-float :complex-double-float))
+ (setq ffi-var (scat "ADDR-REAL-CAST-" var))
+ (nconsc ref-vars
+ `((,ffi-var ,(second (->cffi-type type)) :count 2 :initial-contents (list (realpart ,var) (imagpart ,var))))))
+ (t
+ (setq ffi-var (scat "REF-" var))
+ (nconsc ref-vars
+ `((,ffi-var ,@(->cffi-type type) :initial-element ,var)))))))
+ ;; Output variables
+ (when (and (output-p style) (not (eq type :string)))
+ (nconsc return-vars
+ `((,ffi-var ,var ,type))))
+ ;; Arguments for the lisp wrapper
+ (unless (eq var hidden-var-name)
+ (nconsc defun-args
+ `(,var)))
+ ;; Arguments for the FFI function
+ (nconsc ffi-args
+ `(,ffi-var))
+ ;; Auxillary arguments for FFI
+ (unless (null aux-var)
+ (nconsc aux-ffi-args
+ `(,aux-var))))))
+ ;;Complex returns through hidden variable.
+ (unless (null hidden-var-name)
+ (nconsc aux-args `((,hidden-var-name ,(ecase (second (first pars))
+ (:complex-single-float #c(0e0 0e0))
+ (:complex-double-float #c(0d0 0d0)))))))
+ ;;Keyword argument list
(unless (null defun-keyword-args)
- (setq defun-keyword-args (append '(&key) defun-keyword-args)))
+ (setq defun-keyword-args (cons '&key defun-keyword-args)))
+ ;;Return the function definition
+ (let ((retvar (gensym)))
+ `(
+ ,(recursive-append
+ `(defun ,name ,(append defun-args defun-keyword-args)
+ ,@doc)
+ ;;
+ (unless (null aux-args)
+ `(let (,@aux-args)))
+ ;;Don't use with-foreign.. if ref-vars is nil
+ (unless (null ref-vars)
+ `(with-foreign-objects-stacked (,@ref-vars)))
+ ;;Don't use with-vector-dat.. if array-vars is nil
+ (unless (null array-vars)
+ `(with-vector-data-addresses (,@array-vars)))
+ ;;Declare callbacks
+ callback-code
+ ;;Call the foreign-function
+ `(let ((,retvar (,ffi-fn ,@ffi-args ,@aux-ffi-args)))
+ ;;Ignore return if type is :void
+ ,@(when (eq return-type :void)
+ `((declare (ignore ,retvar))))
+ ;; Copy values in reference pointers back to local
+ ;; variables. Lisp has local scope; its safe to
+ ;; modify variables in parameter lists.
+ ,@(mapcar #'(lambda (decl)
+ (destructuring-bind (ffi-var var type) decl
+ (if (member type '(:complex-single-float :complex-double-float))
+ `(setq ,var (complex (cffi:mem-aref ,ffi-var ,(second (->cffi-type type)) 0)
+ (cffi:mem-aref ,ffi-var ,(second (->cffi-type type)) 1)))
+ `(setq ,var (cffi:mem-aref ,ffi-var ,@(->cffi-type type))))))
+ (remove-if-not #'(lambda (x)
+ (member (first x) ref-vars :key #'car))
+ return-vars))
+ (values
+ ,@(unless (eq return-type :void)
+ `(,retvar))
+ ,@(mapcar #'second return-vars)))))))))
+
+
+(defun def-fortran-callback (func callback-name return-type parm)
+ (let* ((hack-return-type `,return-type)
+ (hack-parm `(,@parm))
+ (hidden-var-name nil))
+ ;;
+ (when (member hack-return-type '(:complex-single-float :complex-double-float))
+ (setq hidden-var-name (gensym "HIDDEN-COMPLEX-RETURN-"))
+ (setq hack-parm `((,hidden-var-name ,hack-return-type :output)
+ ,@parm))
+ (setq hack-return-type :void))
+ ;;
+ (let* ((new-pars nil)
+ (aux-pars nil)
+ (func-pars nil)
+ (array-vars nil)
+ (return-vars nil)
+ (ref-vars nil))
+ (dolist (decl hack-parm)
+ (destructuring-bind (var type &optional (style :input)) decl
+ (let ((ffi-var nil)
+ (func-var nil))
+ (cond
+ ;; Callbacks are tricky.
+ ((callback-type-p type)
+ (setq ffi-var var)
+ (setq func-var var))
+ ;;
+ ((array-p type)
+ (setq ffi-var (scat "ADDR-" var))
+ (setq func-var var)
+ (nconsc array-vars `((,func-var (make-foreign-vector :pointer ,ffi-var :type ,(second (->cffi-type type))
+ :size ,(if-let (size (get-arg :size type))
+ size
+ 1))))))
+ ;;
+ ((string-p type)
+ (setq ffi-var var)
+ (setq func-var var)
+ (nconsc aux-pars
+ `((,(scat "LEN-" var) ,@(->cffi-type :integer)))))
+ ;;
+ ((eq style :input-value)
+ (setq ffi-var var)
+ (setq func-var var))
+ ;; Pass-by-reference variables
+ (t
+ (cond
+ ((member type '(:complex-single-float :complex-double-float))
+ (setq ffi-var (scat "ADDR-REAL-CAST-" var))
+ (setq func-var var)
+ (nconsc ref-vars
+ `((,func-var (complex (cffi:mem-aref ,ffi-var ,(second (->cffi-type type)) 0)
+ (cffi:mem-aref ,ffi-var ,(second (->cffi-type type)) 1))))))
+ (t
+ (setq ffi-var (scat "REF-" var))
+ (setq func-var var)
+ (nconsc ref-vars
+ `((,func-var (cffi:mem-aref ,ffi-var ,@(->cffi-type type)))))))))
+ ;;
+ (nconsc new-pars `((,ffi-var ,@(get-read-in-type type style))))
+ (nconsc func-pars `(,func-var))
+ (when (and (output-p style) (not (eq type :string)))
+ (nconsc return-vars
+ `((,func-var ,ffi-var ,type)))))))
(let ((retvar (gensym)))
`(
- (defun ,name ,(append defun-args defun-keyword-args)
- ,@doc
- (let (,@(if (not (null hidden-var-name))
- `((,hidden-var-name ,@(if (eq (second (first pars))
- :complex-single-float)
- `(#C(0e0 0e0))
- `(#C(0d0 0d0)))))))
- (with-foreign-objects-stack-ed (,@ref-vars)
- (with-vector-data-addresses (,@array-vars)
- (let* (,@aux-args
- ;;Style warnings are annoying.
- ,@(if (not (eq return-type :void))
- `((,retvar (,ffi-fn ,@ffi-args ,@aux-ffi-args))))
- )
- ,@(if (eq return-type :void)
- `((,ffi-fn ,@ffi-args ,@aux-ffi-args)))
- ;; Copy values in reference pointers back to local
- ;; variables. Lisp has local scope; its safe to
- ;; modify variables in parameter lists.
- ,@(mapcar #'(lambda (decl)
- (destructuring-bind (ffi-var var type) decl
- (if (member type '(:complex-single-float :complex-double-float))
- `(setq ,var (complex (cffi:mem-aref ,ffi-var ,(second (->cffi-type type)) 0)
- (cffi:mem-aref ,ffi-var ,(second (->cffi-type type)) 1)))
- `(setq ,var (cffi:mem-aref ,ffi-var ,@(->cffi-type type))))))
- (remove-if-not #'(lambda (x)
- (member (first x) ref-vars :key #'car))
- return-vars))
- (values
- ,@(if (not (eq return-type :void))
- `(,retvar))
- ,@(mapcar #'second return-vars))))))))))))
+ ,(recursive-append
+ `(cffi:defcallback ,callback-name ,@(get-return-type hack-return-type)
+ (,@new-pars ,@aux-pars))
+ ;;
+ (when ref-vars
+ `(let (,@ref-vars)))
+ ;;
+ (when array-vars
+ `(let (,@array-vars)))
+ ;;
+ `(multiple-value-bind (,retvar ,@(mapcar #'car return-vars)) (funcall ,func ,@func-pars)
+ ,@(when (eq hack-return-type :void)
+ `((declare (ignore ,retvar))))
+ ,@(mapcar #'(lambda (decl)
+ (destructuring-bind (func-var ffi-var type) decl
+ (if (member type '(:complex-single-float :complex-double-float))
+ `(setf (cffi:mem-aref ,ffi-var ,(second (->cffi-type type)) 0) (realpart ,func-var)
+ (cffi:mem-aref ,ffi-var ,(second (->cffi-type type)) 1) (imagpart ,func-var))
+ `(setf (cffi:mem-aref ,ffi-var ,@(->cffi-type type)) ,func-var))))
+ (remove-if-not #'(lambda (x)
+ (member (first x) ref-vars :key #'car))
+ return-vars))
+ ,(if (eq hack-return-type :void)
+ nil
+ retvar))))))))
;; Increment the pointer.
-(defmacro incf-sap (type sap &optional (n 1))
- "Increment the pointer address by one \"slot\"
- depending on the type:
+(defun inc-sap (sap type &optional (n 1))
+"
+ Increment the pointer address by one \"slot\"
+ depending on the type:
:double-float 8 bytes
:single-float 4 bytes
:complex-double-float 8x2 bytes
:complex-single-float 4x2 bytes
- "
- (check-type type symbol)
- `(setf ,sap
- (cffi:inc-pointer ,sap
- ,@(ecase type
- (:double-float `((* ,n 8)))
- (:single-float `((* ,n 4)))
- (:complex-double-float `((* ,n 16)))
- (:complex-single-float `((* ,n 8)))))))
+ "
+ (cffi:inc-pointer sap
+ (ecase type
+ (:double-float (* n 8))
+ (:single-float (* n 4))
+ (:complex-double-float (* n 16))
+ (:complex-single-float (* n 8)))))
+
+(define-modify-macro incf-sap (type &optional (n 1)) inc-sap)
;; Supporting multidimensional arrays is a pain.
;; Only support types that we currently use.
(deftype matlisp-specialized-array ()
- `(or (simple-array (complex double-float) (*))
- (simple-array (complex single-float) (*))
- (simple-array double-float (*))
+ `(or (simple-array double-float (*))
+ ;;
(simple-array single-float (*))
;;
(simple-array (signed-byte 32) *)
diff --git a/src/gemm.lisp b/src/gemm.lisp
index e4d227c..5ca6f21 100644
--- a/src/gemm.lisp
+++ b/src/gemm.lisp
@@ -123,15 +123,15 @@
st-c ldc
:head-a hd-a :head-b hd-b :head-c hd-c))
(progn
- (when (eq job-a :t) (transpose-i! a))
- (when (eq job-b :t) (transpose-i! b))
+ (when (eq job-a :t) (transpose! a))
+ (when (eq job-b :t) (transpose! b))
;;
(symbol-macrolet
((loop-col
(mlet* ((cs-b (col-stride b) :type fixnum)
(cs-c (col-stride c) :type fixnum)
- (col-b (col! b 0) :type ,matrix-type)
- (col-c (col! c 0) :type ,matrix-type))
+ (col-b (col~ b 0) :type ,matrix-type)
+ (col-c (col~ c 0) :type ,matrix-type))
(dotimes (j nc-c)
(when (> j 0)
(setf (head col-b) (+ (head col-b) cs-b))
@@ -140,8 +140,8 @@
(loop-row
(mlet* ((rs-a (row-stride a) :type fixnum)
(rs-c (row-stride c) :type fixnum)
- (row-a (transpose-i! (row! a 0)) :type ,matrix-type)
- (row-c (transpose-i! (row! c 0)) :type ,matrix-type))
+ (row-a (transpose! (row~ a 0)) :type ,matrix-type)
+ (row-c (transpose! (row~ c 0)) :type ,matrix-type))
(dotimes (i nr-c)
(when (> i 0)
(setf (head row-a) (+ (head row-a) rs-a))
@@ -153,8 +153,8 @@
((< nr-c nc-c) loop-row)
(t loop-col)))
;;
- (when (eq job-a :t) (transpose-i! a))
- (when (eq job-b :t) (transpose-i! b))
+ (when (eq job-a :t) (transpose! a))
+ (when (eq job-b :t) (transpose! b))
)))
c))
;;;;
@@ -246,9 +246,17 @@
;
(defmethod gemm! ((alpha number) (a real-matrix) (b real-matrix)
+ (beta cl:real) (c complex-matrix)
+ &optional (job :nn))
+ (let ((r-c (mrealpart~ c)))
+ (declare (type real-matrix c))
+ (gemm! alpha a b 0d0 r-c job))
+ c)
+
+(defmethod gemm! ((alpha number) (a real-matrix) (b real-matrix)
(beta complex) (c complex-matrix)
&optional (job :nn))
- (let ((r-c (mrealpart c))
+ (let ((r-c (mrealpart~ c))
(c-be (complex-coerce beta)))
(declare (type real-matrix c)
(type complex-double-float c-al))
@@ -266,10 +274,10 @@
(defmethod gemm! ((alpha cl:real) (a real-matrix) (b complex-matrix)
(beta cl:real) (c complex-matrix)
&optional (job :nn))
- (let ((r-b (mrealpart b))
- (i-b (mimagpart b))
- (r-c (mrealpart c))
- (i-c (mimagpart c))
+ (let ((r-b (mrealpart~ b))
+ (i-b (mimagpart~ b))
+ (r-c (mrealpart~ c))
+ (i-c (mimagpart~ c))
(r-al (coerce alpha 'double-float))
(r-be (coerce beta 'double-float)))
(declare (type real-matrix r-b i-b r-c i-c)
@@ -280,10 +288,10 @@
(defmethod gemm! ((alpha complex) (a real-matrix) (b complex-matrix)
(beta cl:real) (c complex-matrix)
&optional (job :nn))
- (let ((r-b (mrealpart b))
- (i-b (mimagpart b))
- (r-c (mrealpart c))
- (i-c (mimagpart c))
+ (let ((r-b (mrealpart~ b))
+ (i-b (mimagpart~ b))
+ (r-c (mrealpart~ c))
+ (i-c (mimagpart~ c))
(r-al (coerce (realpart alpha) 'double-float))
(i-al (coerce (imagpart alpha) 'double-float))
(r-be (coerce beta 'double-float)))
@@ -306,10 +314,10 @@
(defmethod gemm! ((alpha cl:real) (a complex-matrix) (b real-matrix)
(beta cl:real) (c complex-matrix)
&optional (job :nn))
- (let ((r-a (mrealpart a))
- (i-a (mimagpart a))
- (r-c (mrealpart c))
- (i-c (mimagpart c))
+ (let ((r-a (mrealpart~ a))
+ (i-a (mimagpart~ a))
+ (r-c (mrealpart~ c))
+ (i-c (mimagpart~ c))
(r-al (coerce alpha 'double-float))
(r-be (coerce beta 'double-float)))
(declare (type real-matrix r-a i-a r-c i-c)
@@ -320,10 +328,10 @@
(defmethod gemm! ((alpha complex) (a complex-matrix) (b real-matrix)
(beta cl:real) (c complex-matrix)
&optional (job :nn))
- (let ((r-a (mrealpart a))
- (i-a (mimagpart a))
- (r-c (mrealpart c))
- (i-c (mimagpart c))
+ (let ((r-a (mrealpart~ a))
+ (i-a (mimagpart~ a))
+ (r-c (mrealpart~ c))
+ (i-c (mimagpart~ c))
(r-al (coerce (realpart alpha) 'double-float))
(i-al (coerce (imagpart alpha) 'double-float))
(r-be (coerce beta 'double-float)))
diff --git a/src/gemv.lisp b/src/gemv.lisp
index dbe8ce8..4ce561b 100644
--- a/src/gemv.lisp
+++ b/src/gemv.lisp
@@ -115,7 +115,7 @@
;
(defmethod gemv! ((alpha cl:real) (A real-matrix) (x real-matrix)
(beta complex) (y complex-matrix) &optional (job :n))
- (let ((r-y (mrealpart y)))
+ (let ((r-y (mrealpart~ y)))
(declare (type real-matrix r-y))
;; y <- \beta * y
(scal! (complex-coerce beta) y)
@@ -133,12 +133,12 @@
(beta cl:real) (y complex-matrix) &optional (job :n))
(let ((r-be (coerce beta 'double-float))
(r-al (coerce alpha 'double-float))
- (r-y (mrealpart y)))
+ (r-y (mrealpart~ y)))
(declare (type double-float r-be r-al)
(type real-matrix r-y))
;; y <- \beta * y
(scal! r-be y)
- ;; (mrealpart y) <- (mrealpart y) + \alpha * A o x
+ ;; (mrealpart~ y) <- (mrealpart~ y) + \alpha * A o x
(real-double-gemv!-typed r-al A x 1d0 r-y job))
y)
@@ -147,13 +147,13 @@
(let ((r-al (coerce (realpart alpha) 'double-float))
(i-al (coerce (imagpart alpha) 'double-float))
(r-be (coerce beta 'double-float))
- (r-y (mrealpart y))
- (i-y (mimagpart y)))
+ (r-y (mrealpart~ y))
+ (i-y (mimagpart~ y)))
(declare (type double-float r-al i-al r-be)
(type real-matrix r-y i-y))
- ;; (mrealpart y) <- \beta * (mrealpart y) + (realpart \alpha) . A o x
+ ;; (mrealpart~ y) <- \beta * (mrealpart~ y) + (realpart \alpha) . A o x
(real-double-gemv!-typed r-al A x r-be r-y job)
- ;; (mimagpart y) <- \beta * (mimagpart y) + (imagpart \alpha) . A o x
+ ;; (mimagpart~ y) <- \beta * (mimagpart~ y) + (imagpart \alpha) . A o x
(real-double-gemv!-typed i-al A x r-be i-y job))
y)
@@ -167,26 +167,26 @@
(defmethod gemv! ((alpha cl:real) (A real-matrix) (x complex-matrix)
(beta cl:real) (y complex-matrix) &optional (job :n))
- (let ((r-x (mrealpart x))
- (i-x (mimagpart x))
- (r-y (mrealpart y))
- (i-y (mimagpart y))
+ (let ((r-x (mrealpart~ x))
+ (i-x (mimagpart~ x))
+ (r-y (mrealpart~ y))
+ (i-y (mimagpart~ y))
(r-al (coerce (realpart alpha) 'double-float))
(r-be (coerce beta 'double-float)))
(declare (type double-float r-al r-be)
(type real-matrix r-x i-x r-y i-y))
- ;; (mrealpart y) <- \beta * (mrealpart y) + \alpha . A o (mrealpart x)
+ ;; (mrealpart~ y) <- \beta * (mrealpart~ y) + \alpha . A o (mrealpart~ x)
(real-double-gemv!-typed r-al A r-x r-be r-y job)
- ;; (mimagpart y) <- \beta * (mimagpart y) + \alpha . A o (mrealpart x)
+ ;; (mimagpart~ y) <- \beta * (mimagpart~ y) + \alpha . A o (mrealpart~ x)
(real-double-gemv!-typed r-al A i-x r-be i-y job))
y)
(defmethod gemv! ((alpha complex) (A real-matrix) (x complex-matrix)
(beta cl:real) (y complex-matrix) &optional (job :n))
- (let ((r-x (mrealpart x))
- (i-x (mimagpart x))
- (r-y (mrealpart y))
- (i-y (mimagpart y))
+ (let ((r-x (mrealpart~ x))
+ (i-x (mimagpart~ x))
+ (r-y (mrealpart~ y))
+ (i-y (mimagpart~ y))
(r-al (coerce (realpart alpha) 'double-float))
(i-al (coerce (imagpart alpha) 'double-float))
(r-be (coerce beta 'double-float)))
@@ -209,26 +209,26 @@
(defmethod gemv! ((alpha cl:real) (A complex-matrix) (x real-matrix)
(beta cl:real) (y complex-matrix) &optional (job :n))
- (let ((r-A (mrealpart A))
- (i-A (mimagpart A))
- (r-y (mrealpart y))
- (i-y (mimagpart y))
+ (let ((r-A (mrealpart~ A))
+ (i-A (mimagpart~ A))
+ (r-y (mrealpart~ y))
+ (i-y (mimagpart~ y))
(r-al (coerce (realpart alpha) 'double-float))
(r-be (coerce beta 'double-float)))
(declare (type double-float r-al r-be)
(type real-matrix r-A i-A r-y i-y))
- ;; (mrealpart y) <- \beta * (mrealpart y) + \alpha . A o (mrealpart x)
+ ;; (mrealpart~ y) <- \beta * (mrealpart~ y) + \alpha . A o (mrealpart~ x)
(real-double-gemv!-typed r-al r-A x r-be r-y job)
- ;; (mimagpart y) <- \beta * (mimagpart y) + \alpha . A o (mrealpart x)
+ ;; (mimagpart~ y) <- \beta * (mimagpart~ y) + \alpha . A o (mrealpart~ x)
(real-double-gemv!-typed r-al i-A x r-be i-y job))
y)
(defmethod gemv! ((alpha complex) (A complex-matrix) (x real-matrix)
(beta cl:real) (y complex-matrix) &optional (job :n))
- (let ((r-A (mrealpart A))
- (i-A (mimagpart A))
- (r-y (mrealpart y))
- (i-y (mimagpart y))
+ (let ((r-A (mrealpart~ A))
+ (i-A (mimagpart~ A))
+ (r-y (mrealpart~ y))
+ (i-y (mimagpart~ y))
(r-al (coerce (realpart alpha) 'double-float))
(i-al (coerce (imagpart alpha) 'double-float))
(r-be (coerce beta 'double-float)))
diff --git a/src/norm.lisp b/src/norm.lisp
index 7087dbe..8298da4 100644
--- a/src/norm.lisp
+++ b/src/norm.lisp
@@ -158,7 +158,7 @@
(declare (type fixnum j))
(setq nrm (max nrm
(with-vector-data-addresses ((addr-store store))
- (incf-sap :double-float addr-store (* j n))
+ (incf-sap addr-store :double-float (* j n))
(dasum n addr-store 1)))))
nrm))
((2 :2) (multiple-value-bind (up sigma vp status)
@@ -173,7 +173,7 @@
(declare (type fixnum i))
(setq nrm (max nrm
(with-vector-data-addresses ((addr-store store))
- (incf-sap :double-float addr-store i)
+ (incf-sap addr-store :double-float i)
(dasum m addr-store n)))))
nrm))
((:f :fro :frob :frobenius)
@@ -181,7 +181,7 @@
(dotimes (j m)
(declare (type fixnum j))
(incf nrm (with-vector-data-addresses ((addr-store store))
- (incf-sap :double-float addr-store (* j n))
+ (incf-sap addr-store :double-float (* j n))
(ddot m addr-store 1 addr-store 1))))
(sqrt nrm)))
(t (error "don't know how to take a ~a-norm of a matrix" p))
@@ -228,7 +228,7 @@
(declare (type fixnum j))
(setq nrm (max nrm
(with-vector-data-addresses ((addr-store store))
- (incf-sap :complex-double-float addr-store (* j n))
+ (incf-sap addr-store :complex-double-float (* j n))
(dzasum n addr-store 1)))))
nrm))
((2 :2) (multiple-value-bind (up sigma vp status)
@@ -243,7 +243,7 @@
(declare (type fixnum i))
(setq nrm (max nrm
(with-vector-data-addresses ((addr-store store))
- (incf-sap :complex-double-float addr-store i)
+ (incf-sap addr-store :complex-double-float i)
(dzasum m addr-store n)))))
nrm))
((:f :fro :frob :frobenius)
@@ -251,7 +251,7 @@
(dotimes (j m)
(declare (type fixnum j))
(incf nrm (with-vector-data-addresses ((addr-store store))
- (incf-sap :double-float addr-store (* j n))
+ (incf-sap addr-store :double-float (* j n))
(realpart (zdotc m addr-store 1 addr-store 1)))))
(sqrt nrm)))
(t (error "don't know how to take a ~a-norm of a matrix" p))
diff --git a/src/realimag.lisp b/src/realimag.lisp
index 2b35319..7f59e78 100644
--- a/src/realimag.lisp
+++ b/src/realimag.lisp
@@ -66,12 +66,111 @@
(in-package "MATLISP")
-#+nil (export '(real
- imag))
+(defun mrealpart~ (mat)
+"
+ Syntax
+ ======
+ (MREALPART~ matrix)
+
+ Purpose
+ =======
+ Returns a new SUB-REAL-MATRIX which is the real part of \"matrix\".
+
+ Store is shared with \"matrix\".
+
+ If \"matrix\" is a scalar, returns its real part.
+
+ See IMAG, REALPART, IMAGPART
+"
+
+ (typecase mat
+ (real-matrix mat)
+ (complex-matrix (make-instance 'sub-real-matrix
+ :parent mat :store (store mat)
+ :nrows (nrows mat) :ncols (ncols mat)
+ :row-stride (* 2 (row-stride mat)) :col-stride (* 2 (col-stride mat))
+ :head (* 2 (head mat))))
+ (number (cl:realpart mat))))
+
+(defun mrealpart (mat)
+"
+ Syntax
+ ======
+ (MREALPART matrix)
+
+ Purpose
+ =======
+ Returns a copy of the real part of \"matrix\".
+
+ If \"matrix\" is a scalar, returns its real part.
+
+ See IMAG, REALPART, IMAGPART
+"
+ (typecase mat
+ (real-matrix (copy mat))
+ (complex-matrix (copy (make-instance 'sub-real-matrix
+ :parent mat :store (store mat)
+ :nrows (nrows mat) :ncols (ncols mat)
+ :row-stride (* 2 (row-stride mat)) :col-stride (* 2 (col-stride mat))
+ :head (* 2 (head mat)))))
+ (number (cl:realpart mat))))
+
+(defun mimagpart~ (mat)
+"
+ Syntax
+ ======
+ (MIMAGPART~ matrix)
+
+ Purpose
+ =======
+ Returns a new SUB-REAL-MATRIX which is the imaginary part of \"matrix\".
-(defgeneric real (matrix)
- (:documentation
- "
+ Store is shared with \"matrix\".
+
+ If \"matrix\" is a real-matrix, returns nil.
+
+ If \"matrix\" is a scalar, returns its imaginary part.
+
+ See IMAG, REALPART, IMAGPART
+"
+ (typecase mat
+ (real-matrix nil)
+ (complex-matrix (make-instance 'sub-real-matrix
+ :parent mat :store (store mat)
+ :nrows (nrows mat) :ncols (ncols mat)
+ :row-stride (* 2 (row-stride mat)) :col-stride (* 2 (col-stride mat))
+ :head (+ 1 (* 2 (head mat)))))
+ (number (cl:imagpart mat))))
+
+
+(defun mimagpart (mat)
+"
+ Syntax
+ ======
+ (MIMAGPART~ matrix)
+
+ Purpose
+ =======
+ Returns a copy of the imaginary part of \"matrix\".
+
+ If \"matrix\" is a scalar, returns its imaginary part.
+
+ See IMAG, REALPART, IMAGPART
+"
+
+ (typecase mat
+ (real-matrix (make-real-matrix-dim (nrows mat) (ncols mat)))
+ (complex-matrix (copy (make-instance 'sub-real-matrix
+ :parent mat :store (store mat)
+ :nrows (nrows mat) :ncols (ncols mat)
+ :row-stride (* 2 (row-stride mat)) :col-stride (* 2 (col-stride mat))
+ :head (+ 1 (* 2 (head mat))))))
+ (number (cl:imagpart mat))))
+
+
+(declaim (inline real))
+(defun real (matrix)
+"
Syntax
======
(REAL matrix)
@@ -82,11 +181,12 @@
If MATRIX is a scalar, returns its real part.
See IMAG, REALPART, IMAGPART
-"))
+"
+ (mrealpart matrix))
+
-(defgeneric imag (matrix)
- (:documentation
- "
+(defun imag (matrix)
+"
Syntax
======
(IMAG matrix)
@@ -97,71 +197,5 @@
If MATRIX is a scalar, returns its imaginary part.
See REAL, REALPART, IMAGPART
-"))
-
-(defmethod real ((x number))
- (realpart x))
-
-(defmethod real ((mat real-matrix))
- (copy mat))
-
-(defmethod real ((mat complex-matrix))
- (let* ((n (nrows mat))
- (m (ncols mat))
- (nxm (number-of-elements mat))
- (store (store mat))
- (new-store (allocate-real-store nxm)))
- (declare (type fixnum n m nxm)
- (type (complex-matrix-store-type *) store)
- (type (real-matrix-store-type *) new-store))
-
- (dcopy nxm store 2 new-store 1)
-
- (make-instance 'real-matrix :nrows n :ncols m :store new-store)))
-
-(defmethod real ((mat standard-matrix))
- (error "don't know how to take the real part of a STANDARD-MATRIX,
-its element types are unknown"))
-
-(defmethod imag ((x number))
- (imagpart x))
-
-(defmethod imag ((mat real-matrix))
- (let ((n (nrows mat))
- (m (ncols mat)))
- (declare (type fixnum n m))
- (make-real-matrix-dim n m)))
-
-;;#+(or :cmu :sbcl)
-(defmethod imag ((mat complex-matrix))
- (let* ((n (nrows mat))
- (m (ncols mat))
- (nxm (number-of-elements mat))
- (store (store mat))
- (new-store (allocate-real-store nxm)))
- (declare (type fixnum n m nxm)
- (type (complex-matrix-store-type *) store)
- (type (real-matrix-store-type *) new-store))
-
- (with-vector-data-addresses ((addr-store store)
- (addr-new-store new-store))
- (incf-sap :double-float addr-store)
- (dcopy nxm addr-store 2 addr-new-store 1))
-
- (make-instance 'real-matrix :nrows n :ncols m :store new-store)))
-
-
-;;#+:allegro
-#+nil
-(defmethod imag ((mat complex-matrix))
- (let* ((n (nrows mat))
- (m (ncols mat))
- (nxm (number-of-elements mat))
- (imag (make-real-matrix-dim n m)))
- (declare (type fixnum n m nxm))
-
- (dotimes (i nxm)
- (declare (type fixnum i))
- (setf (matrix-ref imag i) (imagpart (matrix-ref mat i))))
-
- imag))
+"
+ (mimagpart matrix))
\ No newline at end of file
diff --git a/src/standard-matrix.lisp b/src/standard-matrix.lisp
index c1ef52e..4e0deab 100644
--- a/src/standard-matrix.lisp
+++ b/src/standard-matrix.lisp
@@ -292,175 +292,4 @@ matrix and a number"))
(defmethod make-load-form ((matrix standard-matrix) &optional env)
"MAKE-LOAD-FORM allows us to determine a load time value for
matrices, for example #.(make-matrix ...)"
- (make-load-form-saving-slots matrix :environment env))
-
-;;
-#+nil(defmethod print-object ((matrix standard-matrix) stream)
- (dotimes (i (nrows matrix))
- (dotimes (j (ncols matrix))
- (format stream "~A " (matrix-ref-2d matrix i j)))
- (format stream "~%")))
-
-;;
-(defun transpose! (matrix)
-"
- Syntax
- ======
- (transpose! matrix)
-
- Purpose
- =======
- Exchange row and column strides so that effectively
- the matrix is destructively transposed in place
- (without much effort).
-"
- (cond
- ((typep matrix 'standard-matrix)
- (progn
- (rotatef (nrows matrix) (ncols matrix))
- (rotatef (row-stride matrix) (col-stride matrix))
- matrix))
- ((typep matrix 'number) matrix)
- (t (error "Don't know how to take the transpose of ~A." matrix))))
-
-(defmacro with-transpose! (matlst &rest body)
- `(progn
- ,@(mapcar #'(lambda (mat) `(transpose! ,mat)) matlst)
- ,@body
- ,@(mapcar #'(lambda (mat) `(transpose! ,mat)) matlst)))
-
-;;
-(defgeneric transpose (matrix)
- (:documentation
-"
- Syntax
- ======
- (transpose matrix)
-
- Purpose
- =======
- Create a new matrix object which represents the transpose of the
- the given matrix.
-
- Store is shared with \"matrix\".
-
- Settable
- ========
- (setf (transpose matrix) value)
-
- is basically the same as
-
- (copy! value (transpose matrix))
-"))
-
-(defun (setf transpose) (value matrix)
- (copy! value (transpose matrix)))
-
-(defmethod transpose ((matrix number))
- matrix)
-
-;;
-(defgeneric sub-matrix (matrix origin dim)
- (:documentation
-"
- Syntax
- ======
- (sub-matrix matrix origin dimensions)
-
- Purpose
- =======
- Create a block sub-matrix of \"matrix\" starting at \"origin\"
- of dimension \"dim\", sharing the store.
-
- origin, dim are lists with two elements.
-
- Store is shared with \"matrix\"
-
- Settable
- ========
- (setf (sub-matrix matrix origin dim) value)
-
- is basically the same as
-
- (copy! value (sub-matrix matrix origin dim))
-"))
-
-(defun (setf sub-matrix) (value matrix origin dim)
- (copy! value (sub-matrix matrix origin dim)))
-
-;;
-(defgeneric row (matrix i)
- (:documentation
-"
- Syntax
- ======
- (row matrix i)
-
- Purpose
- =======
- Returns the i'th row of the matrix.
- Store is shared with \"matrix\".
-
- Settable
- ========
- (setf (row matrix i) value)
-
- is basically the same as
-
- (copy! value (row matrix i))
-"))
-
-(defun (setf row) (value matrix i)
- (copy! value (row matrix i)))
-
-;;
-(defgeneric col (matrix j)
- (:documentation
-"
- Syntax
- ======
- (col matrix j)
-
- Purpose
- =======
- Returns the j'th column of the matrix.
- Store is shared with \"matrix\".
-
- Settable
- ========
- (setf (col matrix j) value)
-
- is basically the same as
-
- (copy! value (col matrix j))
-"))
-
-(defun (setf col) (value matrix j)
- (copy! value (col matrix j)))
-
-;;
-(defgeneric diag (matrix &optional d)
- (:documentation
-"
- Syntax
- ======
- (diag matrix &optional (d 0))
-
- Purpose
- =======
- Returns a row-vector representing the d'th diagonal of the matrix.
- [a_{ij} : j - i = d]
-
- Store is shared with \"matrix\".
-
- Settable
- ========
- (setf (diag matrix d) value)
-
- is basically the same as
-
- (copy! value (diag matrix d))
-"))
-
-(defun (setf diag) (value matrix &optional (d 0))
- (copy! value (diag matrix d)))
\ No newline at end of file
+ (make-load-form-saving-slots matrix :environment env))
\ No newline at end of file
diff --git a/src/standard-matrix.lisp b/src/tensor.lisp
similarity index 75%
copy from src/standard-matrix.lisp
copy to src/tensor.lisp
index c1ef52e..5d3819f 100644
--- a/src/standard-matrix.lisp
+++ b/src/tensor.lisp
@@ -1,5 +1,5 @@
;; Definitions of STANDARD-MATRIX
-(in-package :matlisp)
+;;(in-package :matlisp)
;;
(declaim (inline allocate-integer4-store))
@@ -7,6 +7,11 @@
(eval-when (load eval compile)
(deftype integer4-matrix-element-type ()
'(signed-byte 32))
+
+ (deftype index-type ()
+ 'fixnum)
+ (deftype index-array-type (size)
+ '(simple-array index-type (,size)))
)
(defun allocate-integer4-store (size &optional (initial-element 0))
@@ -17,106 +22,50 @@ integer storage. Default INITIAL-ELEMENT = 0."
:initial-element initial-element))
;;
-(declaim (inline store-indexing))
-(defun store-indexing (row col head row-stride col-stride)
- (declare (type (and fixnum (integer 0)) row col head row-stride col-stride))
- (the fixnum (+ head (the fixnum (* row row-stride)) (the fixnum (* col col-stride)))))
-
-(defun blas-copyable-p (matrix)
- (declare (optimize (safety 0) (speed 3))
- (type (or real-matrix complex-matrix) matrix))
- (mlet* ((nr (nrows matrix) :type fixnum)
- (nc (ncols matrix) :type fixnum)
- (rs (row-stride matrix) :type fixnum)
- (cs (col-stride matrix) :type fixnum)
- (ne (number-of-elements matrix) :type fixnum))
- (cond
- ((or (= nc 1) (= cs (* nr rs))) (values t rs ne))
- ((or (= nr 1) (= rs (* nc cs))) (values t cs ne))
- (t (values nil -1 -1)))))
-
-(defun blas-matrix-compatible-p (matrix &optional (op :n))
- (declare (optimize (safety 0) (speed 3))
- (type (or real-matrix complex-matrix) matrix))
- (mlet* (((rs cs) (slot-values matrix '(row-stride col-stride))
- :type (fixnum fixnum)))
- (cond
- ((= cs 1) (values :row-major rs (fortran-nop op)))
- ((= rs 1) (values :col-major cs (fortran-op op)))
- ;;Lets not confound lisp's type declaration.
- (t (values nil -1 "?")))))
-
-(declaim (inline fortran-op))
-(defun fortran-op (op)
- (ecase op (:n "N") (:t "T")))
-
-(declaim (inline fortran-nop))
-(defun fortran-nop (op)
- (ecase op (:t "N") (:n "T")))
-
-(defun fortran-snop (sop)
- (cond
- ((string= sop "N") "T")
- ((string= sop "T") "N")
- (t (error "Unrecognised fortran-op."))))
-
-;;
-(defclass standard-matrix ()
- ((number-of-rows
- :initarg :nrows
- :initform 0
- :accessor nrows
- :type fixnum
- :documentation "Number of rows in the matrix")
- (number-of-cols
- :initarg :ncols
- :initform 0
- :accessor ncols
- :type fixnum
- :documentation "Number of columns in the matrix")
+(defclass standard-tensor ()
+ ((rank
+ :accessor rank
+ :type index-type
+ :documentation "Rank of the matrix: number of arguments for the tensor")
+ (dimensions
+ :accessor dimensions
+ :initarg :dimensions
+ :type (index-array-type *)
+ :documentation "Dimensions of the vector spaces in which the tensor's arguments reside.")
(number-of-elements
- :initform 0
:accessor number-of-elements
:type fixnum
- :documentation "Total number of elements in the matrix (nrows * ncols)")
+ :documentation "Total number of elements in the tensor.")
;;
(head
:initarg :head
:initform 0
:accessor head
- :type fixnum
+ :type index-type
:documentation "Head for the store's accessor.")
- (row-stride
- :initarg :row-stride
- :accessor row-stride
- :type fixnum
- :documentation "Row stride for the store's accessor.")
- (col-stride
- :initarg :col-stride
- :accessor col-stride
- :type fixnum
- :documentation "Column stride for the store's accessor.")
+ (strides
+ :initarg :strides
+ :accessor strides
+ :type (index-array-type *)
+ :documentation "Strides for accesing elements of the tensor.")
(store-size
:accessor stor...
[truncated message content] |
