[Matlisp-users] Matlisp subset in pure CL

[Nicolas N. <Nic...@iw...>]

Hello,

I have written some CL routines for elementary full matrix arithmetic
(vector operations, matrix multiplication and inversion).  I need this in
my PDE toolbox Femlisp (www.femlisp.org) for moderately sized full matrices
appearing as subblocks of sparse matrices.  I first used Matlisp but was
dissatisfied with:

1. The overhead of the FFI call which is present in both CMUCL and ACL (we
   discussed this here).

2. The overhead of each call to BLAS operations

3. Portability.

Recently, I replaced Matlisp by my own Common Lisp code.  It uses matrix
classes which are parametrized on the element type and methods which are
compiled at runtime adapted to the matrix class.  I think this is an
interesting approach and intend to extend this technique also for sparse
matrices in future versions of Femlisp.  I speculate that this development
should bring my application in a speed range comparable with other FEM
toolboxes on unstructured grids.

Example (P4, 2.4GHz), M+!-N adds two NxN-matrices:

MATLISP (2_0beta-2003-10-14, with ATLAS)

M+!-1: 0.27 MFLOPS
M+!-2: 1.25 MFLOPS
M+!-4: 4.99 MFLOPS
M+!-8: 15.25 MFLOPS
M+!-16: 71.39 MFLOPS
M+!-32: 186.41 MFLOPS
M+!-64: 353.20 MFLOPS
M+!-128: 432.96 MFLOPS
M+!-256: 78.03 MFLOPS
M+!-512: 72.94 MFLOPS

FL.MATLISP:
M+!-1: 13.11 MFLOPS
M+!-2: 41.94 MFLOPS
M+!-4: 98.69 MFLOPS
M+!-8: 167.77 MFLOPS
M+!-16: 197.38 MFLOPS
M+!-32: 203.36 MFLOPS
M+!-64: 203.36 MFLOPS
M+!-128: 203.36 MFLOPS
M+!-256: 79.89 MFLOPS
M+!-512: 79.89 MFLOPS

We see that the overhead for the M+!-operation has been largely reduced,
although the peak performance of ATLAS is twice as large.

Now, what is the point of this message?  I do not think that Matlisp
becomes unnecessary because of FL.MATLISP.  In fact, FL.MATLISP comprises
only a small subset of BLAS/LAPACK now, and I do not intend to implement
all of it or to achieve the ATLAS performance for BLAS Level 2 operations
as GEMM! or GESV!.  But FL.MATLISP achieves a lot with comparatively little
code.  And, although FL.MATLISP is not yet perfectly implemented, I think
that there are several things which could be improved in Matlisp in the
direction which I have taken:

1. I like my idea of parametrized matrix classes very much.  I can declare
   matrices of uniform element type as e.g. (standard-matrix 'double-float)
   which returns (and maybe generates at runtime) a class named
   |(STANDARD-MATRIX DOUBLE-FLOAT)| or arbitrary others.  Maybe this would
   be interesting also for Matlisp, while keeping the alias to REAL-MATRIX.

2. My BLAS methods are defined on STANDARD-MATRIX.  When called with a
   subclass, e.g. |(STANDARD-MATRIX DOUBLE-FLOAT)|, they compile a method
   adapted to the subclass.  Probably, Matlisp could be implemented in a
   similar way by suitably interfacing to DGEMM, ZGEMM, etc.

It would also be a large advantage if FL.MATLISP and MATLISP could
interoperate.  This is already possible at a lower level, because I have
kept the Fortran indexing style.  But it might be reasonable to strive for
interoperability also on the CLOS level.

So, if someone has any ideas in which direction one could pursue this
project I would like to hear them.  My routines can be found in the
directory "femlisp:src;matlisp" (if you install Femlisp from
www.femlisp.org) in the FL.MATLISP package.  It would also be relatively
easy to extract these routines into a separate project, e.g. CL-MATLISP, if
there is large interest in such a move.

Thank you for any suggestions,

Nicolas.