[lisplab-cvs] r19 - src/matrix
Jørn Inge Vestgården
jivestgarden at common-lisp.net
Tue May 12 19:40:37 UTC 2009
Author: jivestgarden
Date: Tue May 12 15:40:36 2009
New Revision: 19
Log:
Added array specializations for operators and functions
Modified:
src/matrix/level1-array.lisp
src/matrix/level2-array-functions.lisp
Modified: src/matrix/level1-array.lisp
==============================================================================
--- src/matrix/level1-array.lisp (original)
+++ src/matrix/level1-array.lisp Tue May 12 15:40:36 2009
@@ -27,6 +27,15 @@
"True for an array of rank 2"
(= (rank a) 1))
+(defmethod copy ((a array))
+ (if (vector? a)
+ (copy-seq a)
+ (let ((y (make-array (dim a) :element-type (element-type a))))
+ (dotimes (i (size a))
+ (setf (row-major-aref y i)
+ (row-major-aref a i)))
+ y)))
+
(defmethod new ((class (eql 'array)) dim &optional (element-type t) (value 0))
(unless (consp dim) (setf dim (list dim 1)))
(make-array dim
Modified: src/matrix/level2-array-functions.lisp
==============================================================================
--- src/matrix/level2-array-functions.lisp (original)
+++ src/matrix/level2-array-functions.lisp Tue May 12 15:40:36 2009
@@ -1,5 +1,6 @@
;;; Lisplab, level2-array-functions.lisp
-;;; Level2, algbra functions on arrays
+;;; Level2, algbraic functions on arrays
+;;;
;;; TOOD: Make fast methods also for integers.
@@ -20,30 +21,43 @@
;;; 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
(in-package :lisplab)
-
+
+;;; Array specialization for the main algebraic operators and functions.
+;;; It is obviously usefull to have .+, .-, etc. and .sin .cos, etc.
+;;; to operate on arrays. But what about predicates: .=, .<, etc.
+;;; It might be an idea to do specialized these for arrays, but so far
+;;; they are not.
+
+
(defmacro define-array-binary-operator (new old)
+ ;; Creates methods for binary operators specialized on arrys.
+ ;; This is not easy to do in an efficient way, which is one
+ ;; main reason to make specialized matrices (such as blas-real).
+ ;; However I have put in optimizations for double floats. It
+ ;; might be usefull to also do so for complex double-float and some
+ ;; integer types and bytes.
(let ((a (gensym))
(b (gensym))
(c (gensym))
(i (gensym)))
`(progn
- ;; two array
+ ;; two arrays
(defmethod ,new ((,a array) (,b array))
(if (and (eql (element-type ,a) 'double-float)
(subtypep (type-of ,a) 'simple-array)
(eql (element-type ,b) 'double-float)
(subtypep (type-of ,b) 'simple-array))
(let ((,c (copy ,a)))
- (declare ((simple-array double-float) ,b ,c))
+ (declare (type (simple-array double-float) ,b ,c))
(dotimes (,i (min (size ,c) (size ,a)))
(setf (row-major-aref ,c ,i)
(,old (row-major-aref ,c ,i) (row-major-aref ,b ,i))))
,c)
- (let ((,c (create ,a t)))
+ (let ((,c (create ,a 0)))
(dotimes (,i (min (size ,c) (size ,a)))
(setf (vref ,c ,i)
- (,new (vref ,c ,i) (vref ,b ,i))))
+ (,new (vref ,a ,i) (vref ,b ,i))))
,c)))
;; array and number
@@ -53,15 +67,15 @@
(realp ,b))
(let ((,b (coerce ,b 'double-float))
(,c (copy ,a)))
- (declare ((simple-array double-float) ,c))
+ (declare (type (simple-array double-float) ,c))
(dotimes (,i (size ,c))
(setf (row-major-aref ,c ,i)
(,old (row-major-aref ,c ,i) ,b)))
,c)
- (let ((,c (create ,a t)))
+ (let ((,c (create ,a 0)))
(dotimes (,i (size ,c))
(setf (vref ,c ,i)
- (,new (vref ,c ,i) ,b)))
+ (,new (vref ,a ,i) ,b)))
,c)))
;; number and array
@@ -71,15 +85,15 @@
(realp ,a))
(let ((,b (coerce ,a 'double-float))
(,c (copy ,b)))
- (declare ((simple-array double-float) ,c))
+ (declare (type (simple-array double-float) ,c))
(dotimes (,i (size ,c))
(setf (row-major-aref ,c ,i)
(,old ,b (row-major-aref ,c ,i))))
,c)
- (let ((,c (create ,b t)))
+ (let ((,c (create ,b 0)))
(dotimes (,i (size ,c))
(setf (vref ,c ,i)
- (,new ,b (vref ,c ,i))))
+ (,new ,b (vref ,b ,i))))
,c))))))
(define-array-binary-operator .add +)
@@ -88,16 +102,148 @@
(define-array-binary-operator .div /)
(define-array-binary-operator .expt expt)
+(defmacro each-array-element-df-to-df (x form)
+ "Applies a form on each element of an array. The form must
+make real output for real arguments"
+ (let ((i (gensym))
+ (y (gensym)))
+ `(if (and (eql (element-type ,x) 'double-float)
+ (subtypep (type-of ,x) 'simple-array))
+ (let ((,y (copy ,x)))
+ (declare (type (simple-array double-float) ,y))
+ (dotimes (,i (size ,y))
+ (let ((,x (row-major-aref ,y ,i)))
+ (declare (type double-float ,x))
+ (setf (row-major-aref ,y ,i) ,form)))
+ ,y)
+ (let ((,y (create ,x 0)))
+ (dotimes (,i (size ,y))
+ (let ((,x (vref ,x ,i)))
+ (setf (vref ,y ,i) ,form)))
+ ,y))))
+
+(defmacro each-array-element-df-to-complex-df (x form)
+ "Applies a form on each element of an array. The form must
+make complex output for real arguments"
+ (let ((i (gensym))
+ (y (gensym)))
+ `(if (and (eql (element-type ,x) 'double-float)
+ (subtypep (type-of ,x) 'simple-array))
+ (let ((,y (make-array (dim ,x) :element-type '(complex double-float))))
+ (declare (type (simple-array (complex double-float)) ,y))
+ (dotimes (,i (size ,y))
+ (let ((,x (row-major-aref ,x ,i)))
+ (declare (type double-float ,x))
+ (setf (row-major-aref ,y ,i) ,form)))
+ ,y)
+ (let ((,y (create ,x 0))) ; TOOD must make sure to allow complex values
+ (dotimes (,i (size ,y))
+ (let ((,x (vref ,x ,i)))
+ (setf (vref ,y ,i) ,form)))
+ ,y))))
+
+;;; Trignometric functions
+
+(defmethod .sin ((x array))
+ (each-array-element-df-to-df x (.sin x)))
+
+(defmethod .cos ((x array))
+ (each-array-element-df-to-df x (.cos x)))
+
+(defmethod .tan ((x array))
+ (each-array-element-df-to-df x (.tan x)))
+
+;;; Hyperbolic functions
+
+(defmethod .sinh ((x array))
+ (each-array-element-df-to-df x (.sinh x)))
+
+(defmethod .cosh ((x array))
+ (each-array-element-df-to-df x (.cosh x)))
+
+(defmethod .tanh ((x array))
+ (each-array-element-df-to-df x (.tanh x)))
+
+(defmethod .log ((x array) &optional base)
+ (each-array-element-df-to-df x (.log x base)))
+
+(defmethod .exp ((x array))
+ (each-array-element-df-to-df x (.exp x)))
+
+;;; Bessel functions
+
+(defmethod .besj (n (x array))
+ (each-array-element-df-to-df x (.besj n x)))
+
+(defmethod .besy (n (x array))
+ (each-array-element-df-to-df x (.besy n x)))
+
+(defmethod .besi (n (x array))
+ (each-array-element-df-to-df x (.besi n x)))
+
+(defmethod .besk (n (x array))
+ (each-array-element-df-to-df x (.besk n x)))
+;;; Hankel functions. NB! These are complex with real arguments.
+(defmethod .besh1 (n (x array))
+ (each-array-element-df-to-complex-df x (.besh1 n x)))
+(defmethod .besh2 (n (x array))
+ (each-array-element-df-to-complex-df x (.besh2 n x)))
+#|
-#|
+
+(defmacro define-array-binary-bool-operator (new old)
+ (let ((a (gensym))
+ (b (gensym))
+ (i (gensym)))
+ `(progn
+
+ ;; two arrays
+ (defmethod ,new ((,a array) (,b array))
+ (if (and (eql (element-type ,a) 'double-float)
+ (subtypep (type-of ,a) 'simple-array)
+ (eql (element-type ,b) 'double-float)
+ (subtypep (type-of ,b) 'simple-array))
+ (let ()
+ (declare (type (simple-array double-float) ,a ,b))
+ (dotimes (,i (min (size ,a) (size ,b)) t)
+ (unless (,old (row-major-aref ,a ,i)
+ (row-major-aref ,b ,i))
+ (return-from ,new nil))))
+ (dotimes (,i (min (size ,a) (size ,b)) t)
+ (unless (,new (vref ,a ,i)
+ (vref ,b ,i))
+ (return-from ,new nil)))))
+
+ ;; array and number
+ (defmethod ,new ((,a array) (,b number))
+ (if (and (eql (element-type ,a) 'double-float)
+ (subtypep (type-of ,a) 'simple-array)
+ (eql (element-type ,b) 'double-float)
+ (subtypep (type-of ,b) 'simple-array))
+ (let ()
+ (declare (type (simple-array double-float) ,a ,b))
+ (dotimes (,i (min (size ,a) (size ,b)) t)
+ (unless (,old (row-major-aref ,a ,i)
+ (row-major-aref ,b ,i))
+ (return-from ,new nil))))
+ (dotimes (,i (min (size ,a) (size ,b)) t)
+ (unless (,new (vref ,a ,i)
+ (vref ,b ,i))
+ (return-from ,new nil)))))
+
+ ;; number and array
+ (defmethod ,new ((,a number) (,b array))))))
+
+(define-array-binary-bool-operator .< <)
+
#+nil (defun combine-types (a b)
(typecase a
More information about the lisplab-cvs
mailing list