Common Lisp is a wonderful language for explorative computing, but often lacks tools for numerical computing. This project is an attempt at plugging this perceived shortcoming.

Matlisp provides an infix reader for evaluating mathematical expressions. We use ‘/’ as both a binary and unary operator within matlisp’s infix reader, this allows us to ensure semantic uniformity even when dealing with non-commutativity rings.

Solving $A x = b$ is consequently as simple as,

MATLISP> #i(/A * b)

, this operation being generic to A being a scalar/matrix/permutation. It’s also extremely simple to define custom syntax by tweaking “src/reader/infix.lisp”.

For a variety of reasons, Matlisp duplicates general BLAS functionality using Lisp code. These are similar to the naive-loop versions found in the BLAS reference implementation, but are generated using macros, and are consequently quite terse. These macro expand into lean code blocks which are often very competitive (within 10-20%) of the corresponding C code. This allows users greater flexibility when dealing with specific problems. For instance GEMM is essentially generated by the following code-snippet (exercise: why ‘j k i’ ?),

> (defun mm (A B C)
    (einstein-sum #.(tensor 'double-float) (j k i) (ref C i j) (* (ref A i j) (ref B j k))))

Matlisp also switches dynamically between BLAS routines (when available) and native lisp code to avoid FFI overheads. This functionality can be used (with minimal tweaking) without the need for a a BLAS library, or for commutative-ring types which don’t have BLAS routines available.

Most classes/methods in Matlisp are dynamically generated for specialized types at run-time. Generating class specialized for certain field-types is not particularly interesting, and takes place along usual lines.

M> (tensor 'double-float 'dense-tensor)
|<BLAS-MIXIN DENSE-TENSOR: DOUBLE-FLOAT>|
M> (tensor 'double-float 'graph-tensor)
|<GRAPH-TENSOR: DOUBLE-FLOAT>|

Tying method generation with CLOS, however, comes with its own set of complications. Consider this not-so-hypothetical situation: we’re given an abstract class dense-tensor, and we wish to write a method for computing the square of a matrix. This sounds simple enough, and as a first pass, one ends up with something like,

(defgeneric sqmm (x)
  (:method :before ((x tensor))
     (assert (typep x 'tensor-square-matrix) nil "X must be square")))

A first solution can be as simple as,

(defmethod sqmm ((x dense-tensor) &aux (cl-x (class-name (class-of x))))
  (compile-and-eval
    `(defmethod sqmm ((x ,cl-x) &aux (ret (zeros (dimensions x) (class-of x))))
	 (einstein-sum ,cl-x (j k i) (ref ret i j) (* (ref x i j) (ref x j k)))))
  (sqmm x))

This works because when a subclass of dense-tensor is encountered which isn’t specialized, a method is generated and this method is subsequently called. The next time sqmm is called with an instance of this class, closer-mop:compute-applicable-methods-using-classes places this method at the top of the list. This scheme works perfectly if dense-tensor only has descendents at distance 1 and no more (why ?). Fixing this within the confines of CLOS is troublesome, because every method then needs to check if it is ‘special enough’. Cue: Metaobjects to the rescue.

This is solved in Matlisp in three steps,

• Define a subclass tensor-method-generator of closer-mop:funcallable-standard-class (a subclass of a function).
• Define subclasses of closer-mop:specializer, one each for marking generated function arguments, and generator arguments.
• Define custom method dispatch closer-mop:compute-applicable-methods-using-classes for this subclass of generic-functions.

These definitions can be found in ‘src;base;generator.lisp’.

The solution to the question posed is then simply,

(defgeneric sqmm (x)
  (:method :before ((x tensor))
     (assert (typep x 'tensor-square-matrix) nil "X must be square"))
  (:generic-function-class tensor-method-generator))
(defmethod sqmm ((x #.(make-instance 'group-specializer :object-class 'dense-tensor :group-name :x)) &aux (cl-x (class-name (class-of x))))
  (compile-and-eval
    `(defmethod sqmm ((x ,(make-instance 'classp-specializer :object-class (find-class cl-x))) &aux (ret (zeros (dimensions x) ',cl-x)))
	 (einstein-sum ,cl-x (j k i) (ref ret i j) (* (ref x i j) (ref x j k)))))
  (sqmm x))

Methods with classp-specializer are only selected by closer-mop:compute-applicable-methods-using-classes if the argument class match that of the classp-specializer instances which are passed as specializers of the method. Similarly, the group-specializer is only selected if all arguments are subclasses of object-class, and all argument classes for a given group-name are the same. The CLOS discriminator also ensures that the latter method has very low precedence. This solves the problem posed earlier and code-generating ‘Templates’ essentially end up becoming first-class functions.

For debuggability and other reasons, this mechanism is abstracted away inside the macro matlisp:define-tensor-method.

• CFFI

  Matlisp uses CFFI for interfacing with foreign libraries. CFFI support varies between lisp implementations; Matlisp assumes that certain simple-vectors in the lisp implementation are in unboxed form, and have a SAP that can be passed to Fortran functions. It might be possible to replace this with static/pinned arrays, as is kind-of done a the moment for foreign dense-tensor.

  Most of the BLAS functionality is duplicated in native lisp code, and so one can very possibly run Matlisp without any foreign libraries (but at the cost of losing all LAPACK functionality).

  CFFI is provided by quicklisp.

• closer-MOP

  Matlisp now uses the Meta-Object protocol to generate and maintain classes & methods. This is not part of the ANSI standard, and consequently only seems to work on SBCL at the moment. Your luck with other implementations may wary.

• Trivia

  You may also want to use the latest version of Trivia pattern matching library. Matlisp uses the λlist matcher which is a recent addition to the project.

Matlisp has been successfully built on SBCL. Getting it running is thematic of quicklisp,

• Install quicklisp http://www.quicklisp.org/beta/.
• Download matlisp source:

> git clone [email protected]:matlisp/matlisp.git
> ln -s $PWD/matlisp <quicklisp-directory>/local-projects

Fire up your lisp implementation and load as usual with quicklisp:

CL-USER> (ql:quickload :matlisp)
CL-USER> (in-package :matlisp)
M> (named-readtables:in-readtable :infix-dispatch-table)

Matlisp searches for “libblas” (and possibly “liblapack”) first in the directory where the ASDF definition file resides, and then in /usr/lib/ (on linux) and then in paths defined by ld.conf. Linking to something like MKL is as simple as placing a symbolic link by name “libblas.so” in the top-level directory. See src/package.lisp for details. Matlisp relies only on the Fortran API (not CBLAS).

;;Creation
M> (copy! (randn '(2 2)) (zeros '(2 2) '((complex double-float))))
#<|(COMPLEX DOUBLE-FLOAT) STRIDE-ACCESSOR SIMPLE-ARRAY| #(2 2)
  0.8492   -1.976
  2.207    -1.251
>
;;gemv
M> (let ((a (randn '(2 2)))
	   (b (randn 2)))
     #i(a * b))
#<|DOUBLE-FLOAT STRIDE-ACCESSOR SIMPLE-ARRAY| #(2)
1.1885     0.95746
>

;;Tensor contraction
M> (let ((H (randn '(2 2 2)))
	   (b (randn 2))
	   (c (randn 2))
	   (f (zeros 2)))
	   ;;#i(H @ b @ c)
     (einstein-sum #.(tensor 'double-float) (i j k) (ref f i) (* (ref H i j k) (ref b j) (ref c k))))
#<|DOUBLE-FLOAT STRIDE-ACCESSOR SIMPLE-ARRAY| #(2)
0.62586     -1.1128
>
;;Similarly
M> (let ((H (randn '(2 2 2))))
     #i(H @ randn(2) @ randn(2)))
#<|DOUBLE-FLOAT STRIDE-ACCESSOR SIMPLE-ARRAY| #(2)
 0.3234  -0.6201
>

More documentation will be added as the project reaches a stable state.

• matlisp-forbi

  The API for BLAS functions dot ensures inconsistent ABIs between compilers. This package provides a Fortran wrapper (and Lisp methods for `dot`) that fixes these issues. It also provides F77 methods for elementwise division, which follow the `scal` API.e

• Weyl

  Weyl is a CAS written in Lisp (and for Lisp!) at Cornell by Richard Zippel’s group. Currently, this used only within ‘src;base;symbolic.lisp’ (and assoc. infix readers), for working with symbolic expressions. In order to use this functionality, Weyl must be loaded before Matlisp.

  Weyl can installed from ‘[email protected]:matlisp/weyl.git’.

• Dynamic class/method generation, using MOP
• Complete BLAS/LAPACK functionality for types double-float, single-float, (complex single-float), (complex double-float).
• Partial support for dense-tensor with a foreign-pointer store.
• Inplace slicing, real-imaginary part views, negative strides for dense-tensors.
• permutations, sorting, conversion between action/cycle/flip representations.
• Optimized loop generators (einstein/iterate for-mod) in Lisp; BLAS functionality duplicated, and switches automatically b/w Lisp and Fortran.
• Arbitrary tensor contraction.
• Graphs: a general CCS/CCR matrix implementation, lisp adjacency list support, iterate macros for DFS/BFS/SFD graph-traversal, tree-decomposition, cholesky-covers, maximum acyclic subgraph, Djikstra’s algorithms.
• Data structures: Fibonacci heap, Union-Find structure, minimal Doubly linked lists.
• Hash-table sparse tensor: O(1) read/write.
• Co-ordinate sparse tensor

• Unify slicing syntax Unify the slicing syntax used by iterate for-mod/einstein macros. Unify these with a more powerful language.
• Automatic Differentiation
• Symbolic Integration Needs extensive hacking of Weyl.
• Gnuplot interface
• (C)Python-bridge
• FiveAM tests

Matlisp uses a variety of Unicode symbols for some function names and certain operators in the infix reader. The user can readily change these to his suiting, or instead use the following Emacs shortcuts to enter these characters.

;; Lisp
(defun add-lisp-slime-hook (func)
  (add-hook 'lisp-mode-hook func)
  (add-hook 'slime-repl-mode-hook func))
;;#\GREEK_SMALL_LETTER_LAMDA is bound to lambda in :infix-dispatch-table; inherited from λ-reader
(add-lisp-slime-hook #'(lambda () (local-set-key (kbd "C-c \\") (lambda () (interactive (insert "λ"))))))
;;#\LONG_RIGHTWARDS_ARROW_FROM_BAR used for anonymous function definition in Infix
;;#i([x] ⟼ x + 1)
(add-lisp-slime-hook #'(lambda () (local-set-key (kbd "C-c /") (lambda () (interactive (insert "⟼"))))))
;;#\CIRCLE_TIMES used for tensor-product in Infix
(add-lisp-slime-hook #'(lambda () (local-set-key (kbd "C-c *") (lambda () (interactive (insert "⊗"))))))
;;#\MIDDLE_DOT used for tensor-contraction (also bound is @) in Infix
(add-lisp-slime-hook #'(lambda () (local-set-key (kbd "C-c .") (lambda () (interactive (insert "·"))))))
;;Used in the function `(δ-i g &optional i j)` in graph-accessor.lisp
(add-lisp-slime-hook #'(lambda () (local-set-key (kbd "C-c a") (lambda () (interactive (insert "δ"))))))
;;#\DEVANAGARI_LETTER_SA used in infix for tensors involving symbolic expressions.
(add-lisp-slime-hook #'(lambda () (local-set-key (kbd "C-c s") (lambda () (interactive (insert "स"))))))