[lisplab-cvs] r117 - src/matrix
Jørn Inge Vestgården
jivestgarden at common-lisp.net
Thu Dec 10 20:18:33 UTC 2009
Author: jivestgarden
Date: Thu Dec 10 15:18:32 2009
New Revision: 117
Log:
optimizations for cross real and complex. Untested
Modified:
src/matrix/level2-matrix-zge.lisp
src/matrix/store-operators.lisp
Modified: src/matrix/level2-matrix-zge.lisp
==============================================================================
--- src/matrix/level2-matrix-zge.lisp (original)
+++ src/matrix/level2-matrix-zge.lisp Thu Dec 10 15:18:32 2009
@@ -21,6 +21,7 @@
(let ((rx (coerce (realpart value) 'double-float))
(cx (coerce (imagpart value) 'double-float))
(store (matrix-store a)))
+ (declare (type type-blas-store store))
(loop for i from 0 below (length store) by 2 do
(setf (aref store i) rx
(aref store (1+ i)) cx))))
@@ -92,7 +93,7 @@
(expand-generic-function-cdfa-cdf-map)
-;;; Real and matrix
+;;; Complex and matrix
(define-constant +generic-function-cdf-cdfa-map+
'((.add . +_cdf-cdfa)
@@ -167,8 +168,8 @@
`(progn
(def-cross-complex-real-method ,name complex matrix-base-dge)
(def-cross-complex-real-method ,name matrix-base-dge complex)
- (def-cross-complex-real-method ,name matrix-base-zge matrix-base-dge)
- (def-cross-complex-real-method ,name matrix-base-dge matrix-base-zge)
+ ;; (def-cross-complex-real-method ,name matrix-base-zge matrix-base-dge)
+ ;; (def-cross-complex-real-method ,name matrix-base-dge matrix-base-zge)
'done))
(def-all-cross-complex-real-methods .add)
@@ -177,6 +178,82 @@
(def-all-cross-complex-real-methods .div)
(def-all-cross-complex-real-methods .expt)
+;;; Add
+
+(defmethod .add ((a matrix-base-zge) (b real))
+ (let ((c (mcreate a)))
+ (+_cdfa-df (matrix-store a) (coerce b 'double-float) (matrix-store c))
+ c))
+
+(defmethod .add ((b real) (a matrix-base-zge) )
+ (.add a b))
+
+(defmethod .add ((a matrix-base-zge) (b matrix-base-dge))
+ (let ((c (mcreate a)))
+ (+_cdfa-dfa (matrix-store a) (matrix-store b) (matrix-store c))
+ c))
+
+(defmethod .add ((b matrix-base-dge) (a matrix-base-zge))
+ (.add a b))
+
+;;; Sub
+
+(defmethod .sub ((a matrix-base-zge) (b real))
+ (let ((c (mcreate a)))
+ (-_cdfa-df (matrix-store a) (coerce b 'double-float) (matrix-store c))
+ c))
+
+(defmethod .sub ((a real) (b matrix-base-zge))
+ (let ((c (mcreate b)))
+ (-_df-cdfa (coerce a 'double-float) (matrix-store b) (matrix-store c))
+ c))
+
+(defmethod .sub ((a matrix-base-zge) (b matrix-base-dge))
+ (let ((c (mcreate a)))
+ (-_cdfa-dfa (matrix-store a) (matrix-store b) (matrix-store c))
+ c))
+
+(defmethod .sub ((a matrix-base-dge) (b matrix-base-zge))
+ (let ((c (mcreate b)))
+ (-_dfa-cdfa (matrix-store a) (matrix-store b) (matrix-store c))
+ c))
+
+;;; Mul
+
+(defmethod .mul ((a matrix-base-zge) (b real))
+ (let ((c (mcreate a)))
+ (*_cdfa-df (matrix-store a) (coerce b 'double-float) (matrix-store c))
+ c))
+
+(defmethod .mul ((b real) (a matrix-base-zge) )
+ (.mul a b))
+
+(defmethod .mul ((a matrix-base-zge) (b matrix-base-dge))
+ (let ((c (mcreate a)))
+ (*_cdfa-dfa (matrix-store a) (matrix-store b) (matrix-store c))
+ c))
+
+(defmethod .mul ((b matrix-base-dge) (a matrix-base-zge))
+ (.mul a b))
+
+;;; Div
+
+(defmethod .div ((a matrix-base-zge) (b real))
+ (let ((c (mcreate a)))
+ (/_cdfa-df (matrix-store a) (coerce b 'double-float) (matrix-store c))
+ c))
+
+(defmethod .div ((a matrix-base-zge) (b matrix-base-dge))
+ (let ((c (mcreate a)))
+ (/_cdfa-dfa (matrix-store a) (matrix-store b) (matrix-store c))
+ c))
+
+(def-cross-complex-real-method .div matrix-base-dge matrix-base-zge)
+
+;;; Expt
+
+(def-cross-complex-real-method .expt matrix-base-zge matrix-base-dge)
+(def-cross-complex-real-method .expt matrix-base-dge matrix-base-zge)
;;;; Ordinary functions
@@ -233,200 +310,3 @@
-#|
-
-(defmacro each-element-function-matrix-base-zge (x form)
- "Applies a form on each element of an matrix-base-zge."
- (let ((i (gensym))
- (y (gensym)))
- `(let* ((,y (copy ,x)))
- (declare (type matrix-base-zge ,y))
- (dotimes (,i (size ,y))
- (let ((,x (vref ,y ,i)))
- (declare (type (complex double-float) ,x))
- (setf (vref ,y ,i)
- ,form)))
- ,y)))
-
-(defmacro expand-matrix-zge-num-num ()
- (cons 'progn
- (mapcar (lambda (name)
- `(defmethod ,(car name) ((x matrix-base-zge))
- (each-element-function-matrix-base-zge x (,(cdr name) x))))
- +functions-complex-to-complex+)))
-
-(expand-matrix-zge-num-num)
-
-(defmethod .log ((x matrix-base-zge) &optional base)
- (if base
- (each-element-function-matrix-base-zge x (log x base))
- (each-element-function-matrix-base-zge x (log x))))
-
-;;; Bessel functions
-
-(defmethod .besj (n (x matrix-base-zge))
- (each-element-function-matrix-base-zge x (.besj n x)))
-
-(defmethod .besy (n (x matrix-base-zge))
- (each-element-function-matrix-base-zge x (.besy n x)))
-
-(defmethod .besi (n (x matrix-base-zge))
- (each-element-function-matrix-base-zge x (.besi n x)))
-
-(defmethod .besk (n (x matrix-base-zge))
- (each-element-function-matrix-base-zge x (.besk n x)))
-
-(defmethod .besh1 (n (x matrix-base-zge))
- (each-element-function-matrix-base-zge x (.besh1 n x)))
-
-(defmethod .besh2 (n (x matrix-base-zge))
- (each-element-function-matrix-base-zge x (.besh2 n x)))
-
-
-|#
-
-#|
-
-#+nil (defmethod .sqr ((x matrix-base-zge))
- (each-element-function-matrix-base-zge x (* x x)))
-
-
-(defmacro expand-on-matrix-zge-lisplab-two-argument-functions-alist ()
- (cons 'progn
- (mapcar (lambda (name)
- `(def-binary-op-matrix-base-zge ,(car name) ,(cdr name)))
- +lisplab-two-argument-functions-alist+)))
-
-(expand-on-matrix-zge-lisplab-two-argument-functions-alist)
-|#
-#|
-;;; Old code
-
-(defmacro def-binary-op-matrix-base-zge (new old)
- ;;; TODO speed up for real numbers. Is it worth the work?
- (let ((a (gensym "a"))
- (b (gensym "b"))
- (len (gensym "len"))
- (store (gensym "store"))
- (store2 (gensym "store2"))
- (i (gensym "i")))
- `(progn
- (defmethod ,new ((,a matrix-base-zge) (,b number))
- (let* ((,a (copy ,a))
- (,store (matrix-store ,a))
- (,b (coerce ,b '(complex double-float)))
- (,len (size ,a)))
- (declare (type (complex double-float) ,b)
- (type type-blas-store ,store)
- (type type-blas-idx ,len))
- (dotimes (,i ,len)
- (setf (ref-blas-complex-store ,store ,i 0 ,len)
- (,old (ref-blas-complex-store ,store ,i 0 ,len) ,b)))
- ,a))
- (defmethod ,new ((,a number) (,b matrix-base-zge))
- (let* ((,b (copy ,b))
- (,store (matrix-store ,b))
- (,a (coerce ,a '(complex double-float)))
- (,len (size ,b)))
- (declare (type (complex double-float) ,a)
- (type type-blas-store ,store)
- (type type-blas-idx ,len))
- (dotimes (,i ,len)
- (setf (ref-blas-complex-store ,store ,i 0 ,len)
- (,old ,a (ref-blas-complex-store ,store ,i 0 ,len))))
- ,b))
- (defmethod ,new ((,a matrix-base-zge) (,b matrix-base-zge))
- (let* ((,a (copy ,a))
- (,store (matrix-store ,a))
- (,store2 (matrix-store ,b))
- (,len (size ,a)))
- (declare (type type-blas-store ,store)
- (type type-blas-store ,store2)
- (type type-blas-idx ,len))
-
- (dotimes (,i ,len)
- (setf (ref-blas-complex-store ,store ,i 0 ,len)
- (,old (ref-blas-complex-store ,store ,i 0 ,len)
- (ref-blas-complex-store ,store2 ,i 0 ,len))))
- ,a))
- (defmethod ,new ((,a matrix-base-zge) (,b matrix-base-dge))
- (let* ((,a (copy ,a))
- (,store (matrix-store ,a))
- (,store2 (matrix-store ,b))
- (,len (size ,a)))
- (declare (type type-blas-store ,store)
- (type type-blas-store ,store2)
- (type type-blas-idx ,len))
- (dotimes (,i ,len)
- (setf (ref-blas-complex-store ,store ,i 0 ,len)
- (,old (ref-blas-complex-store ,store ,i 0 ,len)
- (aref ,store2 ,i))))
- ,a))
- (defmethod ,new ((,a matrix-base-dge) (,b matrix-base-zge))
- (let* ((,b (copy ,b))
- (,store (matrix-store ,a))
- (,store2 (matrix-store ,b))
- (,len (size ,a)))
- (declare (type type-blas-store ,store)
- (type type-blas-store ,store2)
- (type type-blas-idx ,len))
- (dotimes (,i ,len)
- (setf (ref-blas-complex-store ,store2 ,i 0 ,len)
- (,old (aref ,store ,i)
- (ref-blas-complex-store ,store2 ,i 0 ,len))))
- ,b)))))
-
-(def-binary-op-matrix-base-zge .add +)
-
-(def-binary-op-matrix-base-zge .sub -)
-
-(def-binary-op-matrix-base-zge .mul *)
-
-(def-binary-op-matrix-base-zge .div /)
-
-(def-binary-op-matrix-base-zge .expt expt)
-|#
-
-;;; Hm, much shared code here. Could make a unifying macro.
-
-#+nil (defmethod .imagpart ((a matrix-base-zge))
- (let* ((description (create-matrix-description a :et :d))
- (b (make-matrix-instance description (dim a) 0))
- (store-a (matrix-store a))
- (store-b (matrix-store b))
- (len (size a)))
- (declare (type type-blas-store store-a store-b)
- (type type-blas-idx len))
- (dotimes (i len)
- (setf (aref store-b i)
- (aref store-a (1+ (* 2 i)))))
- b))
-
-#+nil (defmethod .realpart ((a matrix-base-zge))
- (let* ((description (create-matrix-description a :et :d))
- (b (make-matrix-instance description (dim a) 0))
- (store-a (matrix-store a))
- (store-b (matrix-store b))
- (len (size a)))
- (declare (type type-blas-store store-a store-b)
- (type type-blas-idx len))
- (dotimes (i len)
- (setf (aref store-b i)
- (aref store-a (* 2 i))))
- b))
-
-#+nil (defmethod .abs ((a matrix-base-zge))
- (let* ((description (create-matrix-description a :et :d))
- (b (make-matrix-instance description (dim a) 0))
- (store-a (matrix-store a))
- (store-b (matrix-store b))
- (len (size a)))
- (declare (type type-blas-store store-a store-b)
- (type type-blas-idx len))
- (dotimes (i len)
- (setf (aref store-b i)
- (let ((x (aref store-a (* 2 i)))
- (y (aref store-a (1+ (* 2 i)))))
- (sqrt (+ (* x x) (* y y))))))
- b))
-
Modified: src/matrix/store-operators.lisp
==============================================================================
--- src/matrix/store-operators.lisp (original)
+++ src/matrix/store-operators.lisp Thu Dec 10 15:18:32 2009
@@ -246,4 +246,101 @@
`(defun-cdfa-cdfa ,(cdr x) ,(car x)))
+operators-cdfa-cdfa-map+)))
-(expand-operators-cdfa-cdfa-map)
\ No newline at end of file
+(expand-operators-cdfa-cdfa-map)
+
+
+;;;; Now, some special cases of real and imaginary mixing
+;;; Other cases could be optimized too, but these cases are not so obvious.
+
+(defun +_cdfa-df (a b c)
+ (declare (type type-blas-store a c)
+ (type double-float b))
+ (dotimes (i (floor (the type-blas-idx (size a)) 2))
+ (let* ((2i (truly-the type-blas-idx (* 2 i)))
+ (2i+1 (the type-blas-idx (1+ 2i))))
+ (declare (type type-blas-idx 2i 2i+1))
+ (setf (aref c 2i) (+ (aref a 2i) b)
+ (aref c 2i+1) (aref a 2i+1) )))
+ c)
+
+(defun +_cdfa-dfa (a b c)
+ (declare (type type-blas-store a b c))
+ (dotimes (i (the type-blas-idx (size b)))
+ (let* ((2i (truly-the type-blas-idx (* 2 i)))
+ (2i+1 (the type-blas-idx (1+ 2i))))
+ (declare (type type-blas-idx 2i 2i+1))
+ (setf (aref c 2i) (+ (aref a 2i) (aref b i))
+ (aref c 2i+1) (aref a 2i+1) )))
+ c)
+
+(defun -_cdfa-df (a b c)
+ (declare (type type-blas-store a c)
+ (type double-float b))
+ (dotimes (i (floor (the type-blas-idx (size a)) 2))
+ (let* ((2i (truly-the type-blas-idx (* 2 i)))
+ (2i+1 (the type-blas-idx (1+ 2i))))
+ (declare (type type-blas-idx 2i 2i+1))
+ (setf (aref c 2i) (- (aref a 2i) b)
+ (aref c 2i+1) (aref a 2i+1) )))
+ c)
+
+(defun -_cdfa-dfa (a b c)
+ (declare (type type-blas-store a b c))
+ (dotimes (i (the type-blas-idx (size b)))
+ (let* ((2i (truly-the type-blas-idx (* 2 i)))
+ (2i+1 (the type-blas-idx (1+ 2i))))
+ (declare (type type-blas-idx 2i 2i+1))
+ (setf (aref c 2i) (- (aref a 2i) (aref b i))
+ (aref c 2i+1) (aref a 2i+1) )))
+ c)
+
+(defun -_df-cdfa (a b c)
+ (declare (type type-blas-store b c)
+ (type double-float a))
+ (dotimes (i (floor (the type-blas-idx (size b)) 2))
+ (let* ((2i (truly-the type-blas-idx (* 2 i)))
+ (2i+1 (the type-blas-idx (1+ 2i))))
+ (declare (type type-blas-idx 2i 2i+1))
+ (setf (aref c 2i) (- a (aref b 2i))
+ (aref c 2i+1) (- (aref b 2i+1) ))))
+ c)
+
+(defun -_dfa-cdfa (a b c)
+ (declare (type type-blas-store a b c))
+ (dotimes (i (the type-blas-idx (size a)))
+ (let* ((2i (truly-the type-blas-idx (* 2 i)))
+ (2i+1 (the type-blas-idx (1+ 2i))))
+ (declare (type type-blas-idx 2i 2i+1))
+ (setf (aref c 2i) (- (aref a i) (aref b 2i))
+ (aref c 2i+1) (- (aref b 2i+1) ))))
+ c)
+
+(defun *_cdfa-df (a b c)
+ ;; Well, the same as +_dfa-df, but length is twice
+ (*_dfa-df a b c))
+
+(defun *_cdfa-dfa (a b c)
+ (declare (type type-blas-store a b c))
+ (dotimes (i (the type-blas-idx (size b)))
+ (let* ((2i (truly-the type-blas-idx (* 2 i)))
+ (2i+1 (the type-blas-idx (1+ 2i))))
+ (declare (type type-blas-idx 2i 2i+1))
+ (setf (aref c 2i) (* (aref a 2i) (aref b i))
+ (aref c 2i+1) (* (aref a 2i+1) (aref b i)))))
+ c)
+
+
+(defun /_cdfa-df (a b c)
+ ;; Well, the same as +_dfa-df, but length is twice
+ (/_dfa-df a b c))
+
+(defun /_cdfa-dfa (a b c)
+ (declare (type type-blas-store a b c))
+ (dotimes (i (the type-blas-idx (size b)))
+ (let* ((2i (truly-the type-blas-idx (* 2 i)))
+ (2i+1 (the type-blas-idx (1+ 2i))))
+ (declare (type type-blas-idx 2i 2i+1))
+ (setf (aref c 2i) (/ (aref a 2i) (aref b i))
+ (aref c 2i+1) (/ (aref a 2i+1) (aref b i)))))
+ c)
+
More information about the lisplab-cvs
mailing list