Hi,
I was cleaning up code that interfaces to BLAS/LAPACK/UMFPack
libraries, and it occured to me that I could unify the interface for
various special matrices using CLOS. The idea is that I create a
matrix class, and specialize that as needed, including dense matrices
which are really just wrappers around 2d arrays.
If possible, I would like dense-matrix element access to be (nearly)
as fast as aref, but I don't know how. I thought of compiler macros,
but I can't think of an advantageous transformation. As a general
issue, I am wondering if it is possible to make "new" datatypes as
fast as "native" ones.
I include some code snippets and benchmarks below, help is
appreciated, I am using SBCL but I am also interested in general
techniques - when I package this up as a library (eventually), I would
prefer if it was general to all implementations.
Thanks,
Tamas
(in-package :cl-user)
;;;; general matrix, will specialize to dense, sparse, lower/upper
;;;; triangular, symmetric, etc
(declaim (optimize (speed 3)))
(defclass matrix ()
((nrow :accessor nrow :initarg :nrow)
(ncol :accessor ncol :initarg :ncol)
(element-type :accessor element-type :initarg :element-type
:initform t)))
(defgeneric mref (matrix i j)
;;; !! make this setf'able
(:documentation "return element [i,j]"))
;;;; dense matrix, simply a wrapper around a 2D array
(defclass dense-matrix (matrix)
((elements :accessor elements)))
(defmethod initialize-instance :after ((matrix dense-matrix) &key)
(with-slots (elements nrow ncol element-type) matrix
(setf elements (make-array (list nrow ncol) :element-type element-type))))
(defmethod mref ((matrix dense-matrix) i j)
(aref (elements matrix) i j))
(defparameter *nrow* 1000)
(defparameter *ncol* 1000)
(defparameter *m* (make-instance 'dense-matrix :nrow *nrow* :ncol *ncol*
:element-type 'double-float))
(defparameter *a* (make-array (list *nrow* *ncol*) :element-type 'double-float))
(time
(let ((sum 0d0))
(dotimes (i *nrow*)
(dotimes (j *ncol*)
(declare (fixnum i j) (double-float sum))
(incf sum (the double-float (mref *m* i j)))))
sum)) ; 0.3s
(time
(let ((sum 0d0))
(declare ((simple-array double-float (* *)) *a*))
(dotimes (i *nrow*)
(dotimes (j *ncol*)
(declare (fixnum i j))
(incf sum (aref *a* i j))))
sum)) ; 0.1s
On Apr 12, 9:44 am, Tamas K Papp <······@gmail.com> wrote:
> I was cleaning up code that interfaces to BLAS/LAPACK/UMFPack
> libraries, and it occured to me that I could unify the interface for
> various special matrices using CLOS. The idea is that I create a
> matrix class, and specialize that as needed, including dense matrices
> which are really just wrappers around 2d arrays.
Tamas,
I've already started a library that does this, but it isn't public.
I'd like to collaborate on this with you, please send me an email and
I will send you my latest documentation. I've sent a copy of the
document to Liam Healy, the maintainer of GSLL, to get his thoughts.
My basic idea is to define the interface, create a reference
implementation entirely in CL and then generate backends to the
foreign functions.
I don't know what you've been using for unit testing, but Liam and I
have been collaborating on floating point extensions to lisp-unit.
They are publicly available at:
http://repo.or.cz/w/lisp-unit.git
Looking forward to working with you on this.
Cheers,
Tom H.
=== Thomas M. Hermann ===
On Apr 12, 6:37 pm, "Thomas M. Hermann" <··········@gmail.com> wrote:
> On Apr 12, 9:44 am, Tamas K Papp <······@gmail.com> wrote:
>
> > I was cleaning up code that interfaces to BLAS/LAPACK/UMFPack
> > libraries, and it occured to me that I could unify the interface for
> > various special matrices using CLOS. The idea is that I create a
> > matrix class, and specialize that as needed, including dense matrices
> > which are really just wrappers around 2d arrays.
>
> Tamas,
>
> I've already started a library that does this, but it isn't public.
> I'd like to collaborate on this with you, please send me an email and
> I will send you my latest documentation. I've sent a copy of the
> document to Liam Healy, the maintainer of GSLL, to get his thoughts.
> My basic idea is to define the interface, create a reference
> implementation entirely in CL and then generate backends to the
> foreign functions.
>
> I don't know what you've been using for unit testing, but Liam and I
> have been collaborating on floating point extensions to lisp-unit.
> They are publicly available at:
>
> http://repo.or.cz/w/lisp-unit.git
>
> Looking forward to working with you on this.
I've got some similar code for lisp-matrix (on github) which starts to
generalize the framework there for symmetric matrices, and there are a
few other specialized types I'll be working on as I need them (have
the statistical modeling framework to perfect, and the types appear as
needed).
Have to admit that I'm using LIFT with some numerical equality
extensions, and am looking at wrapping the GSLL stuff as well (since
for me, it's one of the "big 4", i.e. Tamas' lisp-centric work, rif's
foriegn storage work, and matlisp-gui-alikes being the other 3 that
I'd like to unify).
Of course, I've got my own work, a lisp-centric R replacement with
roots in LispStat, which is driving my numeric needs, so the numerics
are coming slower than I'd like.
On 12 abr, 11:44, Tamas K Papp <······@gmail.com> wrote:
> Hi,
>
> I was cleaning up code that interfaces to BLAS/LAPACK/UMFPack
> libraries, and it occured to me that I could unify the interface for
> various special matrices using CLOS. The idea is that I create a
> matrix class, and specialize that as needed, including dense matrices
> which are really just wrappers around 2d arrays.
>
> If possible, I would like dense-matrix element access to be (nearly)
> as fast as aref, but I don't know how. I thought of compiler macros,
> but I can't think of an advantageous transformation. As a general
> issue, I am wondering if it is possible to make "new" datatypes as
> fast as "native" ones.
I don't know if this is a great idea and to what extent this would
work, but you can do something more or less the way CFFI does it:
(define-compiler-macro mref (&whole whole matrix i j)
(or (mref-expand (get-declared-type-of-or-give-up matrix) matrix i
j)
whole))
(defgeneric mref-expand (matrix-type matrix i j)
(:method ((matrix-type t) matrix i j)
nil))
(defmethod mref-expand ((matrix-type (eql 'dense-matrix)) matrix i j)
`(aref (elements ,matrix) ,i ,j))
The problem would be to get the declared type of some element. There
is http://common-lisp.net/project/parse-declarations/ , but some outer
macro will need to get these declarations and pass them along.
An alternative is to use some implementation-dependant feature (SBCL's
deftransform seems to fit very well for this purpose, just look at a
couple of them from SBCL's source code and you would see how they
work).
>
> I include some code snippets and benchmarks below, help is
> appreciated, I am using SBCL but I am also interested in general
> techniques - when I package this up as a library (eventually), I would
> prefer if it was general to all implementations.
>
> Thanks,
>
> Tamas
>
On Apr 12, 10:44 am, Tamas K Papp <······@gmail.com> wrote:
[...]
> If possible, I would like dense-matrix element access to be (nearly)
> as fast as aref, but I don't know how.
While it's always nice to have this, how necessary is this for the
problems you want to solve with these classes? IME working with
matrices, sparse or dense, it's fairly unusual to do a lot of
accessing of individual elements; I usually want to do something along
the lines of iterating over a rows or columns all at once. I think
you'd be better off from a performance standpoint if you had methods
that mapped functions over rows or columns of matrices (or strides, et
c.), so you could keep as much of the method dispatching as possible
out of your loops.
You would definitely still want to have MREF when you need it, but
your MATRIX-MULTIPLY generic function wouldn't take anything like the
same sort of hit.
Cheers,
Pillsy
On Mon, 13 Apr 2009 09:46:15 -0700, Pillsy wrote:
> On Apr 12, 10:44 am, Tamas K Papp <······@gmail.com> wrote: [...]
>> If possible, I would like dense-matrix element access to be (nearly) as
>> fast as aref, but I don't know how.
>
> While it's always nice to have this, how necessary is this for the
> problems you want to solve with these classes? IME working with
> matrices, sparse or dense, it's fairly unusual to do a lot of accessing
> of individual elements; I usually want to do something along the lines
> of iterating over a rows or columns all at once. I think you'd be better
> off from a performance standpoint if you had methods that mapped
> functions over rows or columns of matrices (or strides, et c.), so you
> could keep as much of the method dispatching as possible out of your
> loops.
>
> You would definitely still want to have MREF when you need it, but your
> MATRIX-MULTIPLY generic function wouldn't take anything like the same
> sort of hit.
That's a good point. Even when I am filling up matrices elementwise
repeatedly (eg value function iteration in economics), the element
access should not be the bottleneck. Thanks for reminding me.
Tamas
Tamas K Papp <······@gmail.com> writes:
> Hi,
>
> I was cleaning up code that interfaces to BLAS/LAPACK/UMFPack
> libraries, and it occured to me that I could unify the interface for
> various special matrices using CLOS. The idea is that I create a
> matrix class, and specialize that as needed, including dense matrices
> which are really just wrappers around 2d arrays.
>
> If possible, I would like dense-matrix element access to be (nearly)
> as fast as aref, but I don't know how. I thought of compiler macros,
> but I can't think of an advantageous transformation. As a general
> issue, I am wondering if it is possible to make "new" datatypes as
> fast as "native" ones.
Perhaps you could use DEFSTRUCT instead.
In the snipped code below, you don't use multiple inheritance, so you
should be able to just use structures instead of classes. You will
still be paying overhead for generic function dispatch, but with
appropriate declarations, you could perhaps get reasonable access to the
underlying array datastructure, since defstruct accessors are generally
a lot faster than CLOS slot accesses.
--
Thomas A. Russ, USC/Information Sciences Institute