[lisplab-cvs] r171 - trunk/src/matrix

[Jørn Inge Vestgården]

Author: jivestgarden
Date: Thu May 20 15:47:06 2010
New Revision: 171

Log:
optimized and fixed .expt

Modified:
   trunk/src/matrix/level2-matrix-dge.lisp
   trunk/src/matrix/store-operators.lisp

Modified: trunk/src/matrix/level2-matrix-dge.lisp
==============================================================================
--- trunk/src/matrix/level2-matrix-dge.lisp	(original)
+++ trunk/src/matrix/level2-matrix-dge.lisp	Thu May 20 15:47:06 2010
@@ -148,7 +148,7 @@
       (.sub . -_dfa-df)
       (.mul . *_dfa-df)
       (.div . /_dfa-df)
-      (.expt . ^_dfa-df) 
+      ;; (.expt . ^_dfa-df) 
       (.max . max_dfa-df) 
       (.min . min_dfa-df)))
 
@@ -176,7 +176,7 @@
       (.sub . -_df-dfa)
       (.mul . *_df-dfa)
       (.div . /_df-dfa)
-      (.expt . ^_df-dfa) 
+      ;; (.expt . ^_df-dfa) 
       (.max . max_df-dfa) 
       (.min . min_df-dfa)))
 
@@ -204,7 +204,7 @@
       (.sub . -_dfa-dfa)
       (.mul . *_dfa-dfa)
       (.div . /_dfa-dfa)
-      (.expt . ^_dfa-dfa) 
+      ;; (.expt . ^_dfa-dfa) 
       (.max . max_dfa-dfa) 
       (.min . min_dfa-dfa)))
 
@@ -225,6 +225,85 @@
 
 (expand-generic-function-dfa-dfa-map)
 
+;;; The expt is an exception, since negative input can give complex exponent
+;;; and since the general case is slow.
+
+(defun all-integer-elements-p (a)
+  "Helper function for .expt"
+  (declare (type-blas-store a))
+  (dotimes (i (length a))
+    (multiple-value-bind (div mod) (ftruncate (aref a i))
+      (declare (ignore div))
+      (unless (zerop mod)
+	(return-from all-integer-elements-p nil))))
+  t)
+	 
+(defmethod .expt ((a matrix-base-dge) (b matrix-base-dge))
+  (cond ((>= (mmin a) 0d0)
+	 (let ((c (mcreate a)))
+	   (^_dfa>=0-dfa (matrix-store a) (matrix-store b) (matrix-store c))
+	   c))
+	((all-integer-elements-p (matrix-store b))
+	 (let ((c (mcreate a)))
+	   (^_dfa-dfa (matrix-store a) (matrix-store b) (matrix-store c))
+	   c))
+	(t (let ((c (mcreate* a :element-type :z))
+		 (a2 (mcreate* a :element-type :z))
+		 (b2 (mcreate* a :element-type :z)))
+	     ;; There should be some better way to upgrade the element type of a matrix 
+	     ;; while keeping the contents.
+	     (copy-contents a a2)
+	     (copy-contents b b2)
+	     ;; This could in theory be a litle more clever since if all exponents 
+	     ;; are even, the output is still real. However, recognizing even integers 
+	     ;; from float input is risky.
+	     (^_cdfa-cdfa (matrix-store a2) (matrix-store b2) (matrix-store c))
+	     c))))
+	  
+(defmethod .expt ((a matrix-base-dge) (b real))
+  "There is a lot of fuzz going on in here. The reason is because
+the important special cases of exponents -3,-2,-1,0,1,2,3 are a factor 10 faster 
+than the general case on SBCL. Furthermor, output can be complex for non-integer exponent." 
+  (multiple-value-bind (div mod) (truncate b)
+    (if (= 0 mod)
+	(let ((c (mcreate a)))
+	  (case div 
+	    (-3 (^_dfa--3 (matrix-store a) (matrix-store c)))
+	    (-2 (^_dfa--2 (matrix-store a) (matrix-store c)))
+	    (-1 (^_dfa--1 (matrix-store a) (matrix-store c)))
+	    (0  (mfill c 1))
+	    (1  (copy-contents a c))
+	    (2  (^_dfa-+2 (matrix-store a) (matrix-store c)))
+	    (3  (^_dfa-+3 (matrix-store a) (matrix-store c)))
+	    (t  (^_dfa-df (matrix-store a) (coerce div 'double-float) (matrix-store c))))
+	  c)
+	(let ((min (mmin a)))
+	  (if (>= min 0)
+	      (let ((c (mcreate a)))
+		(^_dfa>=0-df (matrix-store a) (coerce b 'double-float) (matrix-store c))
+		c)
+	      (let ((c (mcreate* a :element-type :z))
+		    (a2 (mcreate* a :element-type :z)))
+		(copy-contents a a2)
+		(^_cdfa-cdf (matrix-store a2) 
+			    (coerce b '(complex double-float)) 
+			    (matrix-store c))
+		c))))))
+		
+(defmethod .expt ((a real) (b matrix-base-dge) )
+  (if (>= a 0)
+      (let ((c (mcreate b)))
+	(^_df>=0-dfa (coerce a 'double-float) (matrix-store b) (matrix-store c))
+	c)
+      (if (all-integer-elements-p (matrix-store b))
+	  (let ((c (mcreate b)))	      
+	    (^_df-dfa (coerce a 'double-float)  (matrix-store b) (matrix-store c))
+	    c)
+	  (let ((c (mcreate* b :element-type :z))
+		(b2 (mcreate* b :element-type :z)))
+	    (copy-contents b b2)
+	    (^_cdf-cdfa (coerce a '(complex double-float)) (matrix-store b2) (matrix-store c))
+	    c))))
 
 
 (define-constant +generic-function-dfa-to-dfa-map+ ;really bad name    

Modified: trunk/src/matrix/store-operators.lisp
==============================================================================
--- trunk/src/matrix/store-operators.lisp	(original)
+++ trunk/src/matrix/store-operators.lisp	Thu May 20 15:47:06 2010
@@ -132,6 +132,50 @@
 	   (type double-float a))
   (map-into c (lambda (x) (expt a x)) b))
 
+(defun ^_dfa>=0-dfa (a b c)
+  (declare (type (simple-array (double-float 0d0) (*)) a b c))
+  (map-into c #'expt a b))
+
+(defun ^_df>=0-dfa (a b c)
+  (declare (type type-blas-store b c)
+	   (type (double-float 0d0) a))
+  (map-into c (lambda (x) (expt a x)) b))
+
+(defun ^_dfa>=0-df (a b c)
+  (declare (type (simple-array (double-float 0d0) (*)) a c)
+	   (type double-float b))
+  (map-into c (lambda (x) (expt x b)) a))
+
+(defun ^_dfa-+2 (a c)
+  (declare (type type-blas-store a c))
+  (dotimes (i (length a))
+    (setf (aref c i) (expt (aref a i) 2)))
+  c)
+
+(defun ^_dfa-+3 (a c)
+  (declare (type type-blas-store a c))
+  (dotimes (i (length a))
+    (setf (aref c i) (expt (aref a i) 3)))
+  c)
+
+(defun ^_dfa--1 (a c)
+  (declare (type type-blas-store a c))
+  (dotimes (i (length a))
+    (setf (aref c i) (expt (aref a i) -1)))
+  c)
+
+(defun ^_dfa--2 (a c)
+  (declare (type type-blas-store a c))
+  (dotimes (i (length a))
+    (setf (aref c i) (expt (aref a i) -2)))
+  c)
+
+(defun ^_dfa--3 (a c)
+  (declare (type type-blas-store a c))
+  (dotimes (i (length a))
+    (setf (aref c i) (expt (aref a i) -3)))
+  c)
+
 (defun max_dfa-dfa (a b c)
   (declare (type type-blas-store a b c))
   (map-into c #'max a b))