* benchmark/lisp-vs-c-aref.lisp: same as the one on the page of

	Mark Hoemmen but with added tests for the MREF of lisp-matrix and
	matlisp

	* experimental/element-access.lisp: some tries to see if element
	access can be improved while having matrices defined as CLOS
	classes

ChangeLog

@@ -1,3 +1,13 @@
+2008-05-07  Evan Monroig  <evan.monroig@gmail.com>
+
+	* benchmark/lisp-vs-c-aref.lisp: same as the one on the page of
+	Mark Hoemmen but with added tests for the MREF of lisp-matrix and
+	matlisp
+
+	* experimental/element-access.lisp: some tries to see if element
+	access can be improved while having matrices defined as CLOS
+	classes
+
 2008-05-04  Evan Monroig  <evan.monroig@gmail.com>
 
 	* lapack-utils.lisp (def-lapack-method): add wrapper around method

benchmark/lisp-vs-c-aref.lisp

@@ -0,0 +1,160 @@
+;;; Pull out all the optimization stops.
+(declaim (optimize (safety 0) (debug 0) (speed 3)))
+
+;;; Make sure you've loaded ASDF and your CFFI directory was pushed 
+;;; onto asdf:*central-registry*.
+(asdf:oos 'asdf:load-op :cffi)
+
+;;; We set the array size to 4 million so that the benchmark takes a
+;;; decent amount of time to run and so that the array doesn't fit in
+;;; cache (MAKE-ARRAY with the :INITIAL-ELEMENT keyword causes the 
+;;; array to be written to, meaning that the whole array could be in 
+;;; the cache).  Of course if you have an IBM Power5 with its monster
+;;; 36 MB L3 cache, you might want to change this value :P
+(defconstant +ASIZE+ 4000000)
+(declaim (type fixnum +ASIZE+))
+
+
+(defmacro fill-foreign-array (array-name type count with-what)
+  "Fills the given foreign array ARRAY-NAME with COUNT elements 
+   WITH-WHAT, each of which are of (CFFI) type TYPE.  Macro-
+   expansion fails if the macro doesn't recognize the given type."
+  (cond ((or (eq type :double) (eq type 'double))
+	 `(dotimes (i ,count)
+	    (declare (type fixnum i ,count))
+	    (setf (cffi:mem-aref ,array-name :double i) 
+	      (the double-float ,with-what))))
+	((or (eq type :int) (eq type 'int))
+	 `(dotimes (i ,count)
+	    (declare (type fixnum i ,count))
+	    (setf (cffi:mem-aref ,array-name :int i) 
+	      ,with-what)))
+	(t (error "I don't know how to fill with type ~A" type))))
+
+;;; Here is a more general FILL-FOREIGN-ARRAY.  Note that it doesn't 
+;;; do any type checking on ARRAY-NAME.
+#|
+(defmacro fill-foreign-array (array-name type count with-what)
+  `(dotimes (i ,count)
+     (declare (type fixnum i ,count))
+     (setf (cffi:mem-aref ,array-name ,type i) ,with-what)))
+|#   
+
+(defmacro with-foreign-alloc ((array-name type count &optional (init-elt 0 init-elt-supplied-p)) &body body)
+  "Allocates a foreign array (on the heap) named ARRAY-NAME of 
+   CFFI type TYPE containing COUNT elements.  If INIT-ELT is 
+   supplied, all the elements of ARRAY-NAME are set to INIT-ELT.  
+   Then BODY is executed in an UNWIND-PROTECT, and the foreign 
+   array is deallocated."
+  (if init-elt-supplied-p
+      `(let ((,array-name (cffi:foreign-alloc ,type :count ,count)))
+	 (fill-foreign-array ,array-name ,type ,count ,init-elt)
+	 (unwind-protect
+	     (progn
+	       ,@body)
+	   (cffi:foreign-free ,array-name)))
+    `(let ((,array-name (cffi:foreign-alloc ,type :count ,count))
+	   (unwind-protect
+	       (progn
+		 ,@body)
+	     (cffi:foreign-free ,array-name))))))
+
+
+
+(defun lisp-aref-benchmark ()
+  "Benchmark for a 1-D Lisp array, with AREF."
+  (let ((A (make-array +ASIZE+ 
+		       :element-type 'double-float 
+		       :initial-element 1.0d0))
+	(s 0.0d0))
+    (declare (type double-float s))
+    (declare (type (simple-array double-float (*)) A))
+    (dotimes (i +ASIZE+ s)
+       (declare (type fixnum i))
+      (incf s (aref A i)))))
+
+
+
+(defun foreign-aref-benchmark ()
+  "Benchmark for a 1-D C array, with CFFI:MEM-AREF.  The function does pointer
+   arithmetic every time it is called."
+  (with-foreign-alloc (A :double +ASIZE+ 1.0d0)
+    (let ((s 0.0d0))
+      (declare (type double-float s))
+      (dotimes (i +ASIZE+ s)
+	(declare (type fixnum i))
+	(incf s (cffi:mem-aref A :double i))))))
+
+
+(defun optimized-foreign-aref-benchmark ()
+  "Benchmark for a 1-D C array, with CFFI:MEM-REF instead of CFFI:MEM-AREF.
+   We do the pointer arithmetic by hand.  CFFI makes this awkward since a 
+   pointer may be an object different than an (unsigned-byte 32); you have
+   to call INC-POINTER instead of just adding 8 to it.  Thus, the word
+   \"optimized\" in the name of this function should be taken with a grain
+   of salt."
+  (with-foreign-alloc (A :double +ASIZE+ 1.0d0)
+    (let ((s 0.0d0))
+      (declare (type double-float s))
+      ;; We need to make sure that B is a deep copy of A.
+      ;; We'll be playing with B's address but we need to
+      ;; hold on to A so that WITH-FOREIGN-ALLOC will 
+      ;; deallocate the array correctly.
+      (let ((B (cffi:make-pointer (cffi:pointer-address A))))
+	(dotimes (i +ASIZE+ s)
+	  (declare (type fixnum i))
+	  (incf s (cffi:mem-ref B :double))
+	  ;; inc-pointer might cons (it returns a new pointer).
+	  (setf B (cffi:inc-pointer B 8)))))))
+
+(defun lisp-matrix-vref-benchmark ()
+  (let ((a (lisp-matrix::make-vector +ASIZE+ 'lisp-matrix::double :initial-element 1d0))
+        (s 0d0))
+    (declare (type double-float s)
+             (type lisp-matrix::vector-double a))
+    (dotimes (i +ASIZE+ s)
+      (declare (type fixnum i))
+      (incf s (lisp-matrix::vref a i)))))
+
+(defun lisp-matrix-vref-lambda-benchmark ()
+  (declare (optimize (speed 3) (safety 0) (debug 0)))
+  (let ((a (lisp-matrix::make-vector +ASIZE+ 'lisp-matrix::double :initial-element 1d0))
+        (s 0d0))
+    (declare (type double-float s)
+             (type lisp-matrix::vector-double a))
+    (let ((vref (lisp-matrix::vref-lambda a)))
+      (declare (type function vref))
+      (dotimes (i +ASIZE+ s)
+        (declare (type fixnum i))
+        (incf s (the double-float (funcall vref a i)))))))
+
+
+
+(defun matlisp-mref-benchmark ()
+  (let ((a (matlisp::ones +ASIZE+ 1))
+        (s 0d0))
+    (declare (type double-float s)
+             (type matlisp::real-matrix a))
+    (dotimes (i +ASIZE+ s)
+      (declare (type fixnum i))
+      (incf s (matlisp::mref a i)))))
+
+(defun matlisp-mref-benchmark ()
+  (let ((a (matlisp::ones +ASIZE+ 1))
+        (s 0d0))
+    (declare (type double-float s)
+             (type matlisp::real-matrix a))
+    (dotimes (i +ASIZE+ s)
+      (declare (type fixnum i))
+      (incf s (matlisp::mref a i)))))
+
+(defun matlisp-store-aref-benchmark ()
+  (let* ((a (matlisp:ones +ASIZE+ 1))
+         (s 0d0)
+         (st (matlisp::store a)))
+    (declare (type double-float s)
+             (type matlisp::real-matrix a)
+             (type (simple-array double-float (*)) st))
+    (dotimes (i +ASIZE+ s)
+      (declare (type fixnum i))
+      (incf s (aref st i)))))

experimental/element-access.lisp

@@ -0,0 +1,116 @@
+(in-package :lisp-matrix)
+
+(import 'cl-utilities:once-only)
+
+(defmethod vref-code ((a matrix-double) i)
+  (once-only (a i)
+   `((declare (type fixnum ,i))
+     (fnv-double-ref (data ,a) ,i))))
+
+(defmethod vref-lambda ((a matrix-double))
+  (lambda (a i)
+    (declare (type fixnum i))
+    (fnv-double-ref (data a) i)))
+
+(defmethod vref-lambda ((a vector-double))
+  (declare (optimize (speed 3) (safety 0) (debug 0)))
+  (lambda (a i)
+    (declare (type fixnum i)
+             (type vector-double a))
+    (the double-float (fnv-double-ref (data a) i))))
+
+(defmethod mapc-matrix (function (a vector-double))
+  (declare (optimize (speed 3) (safety 0) (debug 0)))
+  (dotimes (i (nelts a))
+    (declare (type fixnum i))
+    (funcall function (fnv-double-ref (data a) i))))
+
+(defun cube (x) (* x x x))
+
+(define-compiler-macro cube (x)
+  (if (numberp x)
+      (* x x x)
+      `(let ((.x ,x))
+         (* .x .x .x))))
+
+(funcall (compiler-macro-function 'cube)
+         '(cube (incf x))
+         nil)
+
+(define-compiler-macro mref (&whole form a i j)
+  (typecase a
+    (matrix-double
+     `(progn
+        (declare (type fixnum ,i ,j))
+        (the double-float
+          (fnv-double-ref (data ,a)
+                          (flatten-matrix-indices ,a ,i ,j)))))
+    (t
+     form)))
+
+
+#+sbcl
+(defun lexenv-get-type (symbol env)
+  (assert (typep symbol 'symbol))
+  (let ((lambda-var
+         (cdr (assoc symbol (sb-c::lexenv-vars env)))))
+    (and lambda-var
+         (typep lambda-var 'sb-c::lambda-var)
+         (sb-c::lambda-var-type lambda-var))))
+
+(defun test-type (form type env)
+  "helper function for compiler macros - test that FORM is of the
+given type, either directly if it is an object, or if there is a (THE
+...) form, or through lexical bindings in the environment ENV."
+  (typecase form
+    (list (and (eql (first form) 'the)
+               (equal (second form) type)
+               (not (null (third form)))))
+    (atom
+     (or (typep form type)
+         (let ((env-type (lexenv-get-type form env)))
+           (and (subtypep env-type type)
+                (subtypep type env-type)))))
+    (t nil)))
+
+(the double-float 1d0)
+
+(define-compiler-macro vref (&whole form x i &environment env)
+  (cond ((test-type x 'vector-double env)
+         (print
+          (let ((fnv-array (gensym "FNV-ARRAY")))
+           `(let ((,fnv-array (data ,x)))
+              (declare (type fixnum ,i)
+                       (type fnv-double ,fnv-array))
+              (the double-float
+                (fnv-double-ref ,fnv-array ,i))))))
+        (t (print form))))
+
+(let ((a (make-matrix 2 2 'double)))
+  (funcall (compiler-macro-function 'mref)
+           `(mref ,a 0 0)
+           nil))
+
+(let ((a (make-matrix 2 2 'float)))
+  (funcall (compiler-macro-function 'mref)
+           `(mref ,a 0 0)
+           nil))
+
+(let ((a (make-vector 2 'double)))
+  (funcall (compiler-macro-function 'vref)
+           `(vref ,a 0)
+           nil))
+
+(let ((a (make-vector 2 'float)))
+  (funcall (compiler-macro-function 'vref)
+           `(vref ,a 0)
+           nil))
+
+
+(funcall (compiler-macro-function 'vref)
+         `(vref ,(make-vector 2 'double) 0)
+         nil)
+
+(let ((a (make-matrix 2 2 'float)))
+  (mref a 1 1))
+

lapack-methods.lisp

@@ -3,7 +3,7 @@
 ;;; This file contains actual LAPACK methods.  See functions in
 ;;; lapack-utils.lisp for how supporting utility macros and functions.
 ;;;
-;;; Time-stamp: <2008-05-04 13:29:49 Evan Monroig>
+;;; Time-stamp: <2008-05-05 21:28:47 Evan Monroig>
 
 (def-lapack-method gemm (alpha (a !matrix-type) (b !matrix-type)
                                &optional (beta 0d0) c)
@@ -35,3 +35,91 @@
                beta
                (data c)
                (real-nrows c))))
+
+
+;; (make-fnv-int32 2 :initial-value 0)
+;; (make-fnv-int32 (ncols a) :initial-value 0)
+
+(defun fnv-type->lisp-type (fnv-type)
+  (ecase fnv-type
+    (double 'double-float)
+    (float 'single-float)
+    (complex-float '(complex single-float))
+    (complex-double '(complex double-float))))
+
+(defmethod rand (nrows ncols fnv-type)
+  "random matrix"
+  (let ((a (make-matrix nrows ncols fnv-type))
+        (one (coerce 1 (fnv-type->lisp-type fnv-type))))
+    (dotimes (i nrows a)
+      (dotimes (j ncols)
+        (setf (mref a i j)
+              (random one))))))
+
+(defmacro call-with-work ((lwork work type) call)
+  `(let ((work (make-matrix 1 1 ,type))
+         (,lwork -1))
+     ,call
+     (setq ,lwork (floor (mref ,work 0 0)))
+     (setq ,work (make-vector ,lwork ,type))
+     ,call))
+
+(defun check-info (info function-name)
+  (unless (= info 0)
+    (error "~a: error in argument ~d" function-name (- info))))
+
+(def-lapack-method gelsy ((a !matrix-type) (b !matrix-type)
+                          rcond &optional jpvt)
+  ;; FIXME: has LWORK and RWORK for %ZGELSY and %CGELSY
+  (unless jpvt
+    (setq jpvt (make-fnv-int32 (ncols a) :initial-value 0)))
+  (let ((rank (make-fnv-int32 1 :initial-value 0))
+        (info (make-fnv-int32 1 :initial-value 0)))
+    ;; FIXME: B needs to be resized anyway if A has more columns than
+    ;; rows, to allow for enough storage for the result matrix =>
+    ;; quick fix: disallow this
+    (assert (<= (ncols a) (nrows a)))
+    (with-copies ((a (or (not unit-stride-p)
+                         (not zero-offset-p)
+                         transposed-p))
+                  (b (or (not unit-stride-p)
+                         (not zero-offset-p)
+                         transposed-p)))
+        ;; FIXME: the value RANK is not correct because the cffi type
+        ;; :FORTRAN-INT does not define a TRANSLATE-FROM-FOREIGN method?
+        ;; => why not use a standard cffi integer type anyway??
+        ;; (a fix is to make RANK a fnv-int32 with one element
+        (progn
+          (check-info (fnv-int32-ref info 0) "GELSY")
+          (list (if (= (nrows b) (ncols a))
+                    b
+                    (window b :nrows (ncols a)))
+                (fnv-int32-ref rank 0)))
+     (call-with-work (lwork work '!data-type)
+                     (!function (nrows a)
+                                (ncols a)
+                                (ncols b)
+                                (data a)
+                                (real-nrows a)
+                                (data b)
+                                (real-nrows b)
+                                jpvt
+                                rcond
+                                rank
+                                (data work)
+                                lwork
+                                info)))))
+
+#+nil
+(let* ((m 10)
+       (n 10)
+       (a (rand m n 'double))
+       (x (rand n 1 'double))
+       (b (gemm 1d0 a x))
+       (rcond (* (coerce (expt 2 -52) 'double-float)
+                 (max (nrows a) (ncols a))))
+       (orig-a (copy a))
+       (orig-b (copy b))
+       (orig-x (copy x)))
+  (list x (gelsy a b rcond)))
+