Compile-time dispatching for fixed-length vector operations

The WITH-FIXED-TYPES macro

The goal is to have generic vector / matrix operators like +, *, etc., that work auto-magically for fixed-length vectors and matrices. However, we don't want to pay runtime CLOS dispatching costs. What we need is a way to tag the types (and sizes) of the vectors and matrices at macro-expansion time so that we can dispatch (at macro-expansion time) from the "generic" (not CLOS) operator to the right expansion.

For example:

(with-fixed-types
    ((dmat4x4 A B C)
     (dvec4 x w))
  (+ (* A B x) (* C w)))
Note that "generic" * needs to know how to optimize these expressions, which may have different combinations of vectors and matrices. (The simplest thing to do is that if the last element in the argument list is a vector, just apply the matrices to it from right to left.)

The idea of WITH-FIXED-TYPES is this:

Optimizing vector operation expressions

Here's the idea: inside of a WITH-FIXED-TYPES, we store the datatype and dimensions (vectors have dimensions (k 1)) of all the type-known objects in the property lists of their symbols. Then we transform expressions with operators into the macro OPAPPLY, which works like this:

(defmacro OPAPPLY (op &rest args) ...)
  1. If any arg in ARGS is a symbol, run TAG-SYMBOL on it to get a tagged expression. Once all the symbols in ARGS have been tagged, run OPAPPLY on the results.
  2. If any arg is an expression with an operator that we know about, recurse on it.
  3. Otherwise, once all the expressions are tagged, combine them according to the rules of the operator.
We tag a symbol like this: Look up its type in the property list (the type was stored there by the WITH-FIXED-TYPES macro). Then tag it like this:
(tag DATATYPE DIMENSIONS-LIST #( (aref x 0) ... ))
That is, we turn the symbol into a symbolic vector / matrix expression. (This is feasible because the vectors and matrices are small.) Combining symbols would work like this:
(tag double-float (2 1) #( (aref x 0) (aref x 1) )) +
(tag double-float (2 1) #( (aref y 0) (aref y 1) ))
->
(tag double-float (2 1) #( (+ (aref x 0) (aref y 0))
                           (+ (aref x 1) (aref y 1)) ))
When the OPAPPLY macro is done, we get a "tagged" list whose last element is an array expression. If the result of OPAPPLY is fed into a DSETV, we don't need to create any temporaries -- just SETF each element of the left-hand side vector with its corresponding element in the right-hand side expression. For example,
(dsetv x (+ y z)) ->
(progn
  (setf (aref x 0) (+ (aref y 0) (aref z 0)))
  (setf (aref x 1) (+ (aref y 1) (aref z 1)))
  x)
(We should of course include appropriate type declarations so that the compiler can optimize everything.) If we don't have a DSETV, then we really want to create a new vector and return it, so we need to figure out how to "strip the tags" and create a new vector of the appropriate specialized type.

We express "operator rules" using Prolog-style matching. For example, + has the rule that all the arguments match both in their datatype and in their dimensions. The * operator is harder because it could take matrices along with vectors. We want the dimensions to match the pattern (n1, n2), (n2, n3), (n3, n4), ... and we also want all the datatypes to match. The * operator should also do automatic parenthesization using the matrix chain algorithm. (You could just implement the matrix chain algorithm naively with memoization.)

The fixed-size types

Names are like this: dvec4, zmat4x4. The first letter is the datatype (using the LAPACK subroutine naming convention: s for single-float, d for double-float, c for (complex single-float), z for (complex double-float)). Then "vec" for vector or "mat" for matrix. Then the dimension(s). This scheme keeps the name short but also conveys all the information in the type.

We should use Lisp arrays to store the data -- 1-d arrays for vectors and 2-d arrays for matrices (for simplicity in coding the macros). We can and should specialize the types of the arrays as far as possible: e.g. a cvec4 is a

(simple-array (complex single-float) 4)
and a dmat4x4 is a
(simple-array double-float (4 4))

Note that Lisp implementations are not required to implement double-float and single-float using the IEEE 754 floating-point types with corresponding names. In fact a Lisp implementation may choose to implement them both using the same datatype. Furthermore, a Lisp implementation need not unbox the elements in a SIMPLE-ARRAY. For example, CLISP boxes floating-point values in a SIMPLE-ARRAY. As a result, since we want to implement the fixed-size vectors and matrices using Lisp arrays, we should implement the operators entirely in Lisp to avoid copy overhead and potential loss of accuracy (e.g. if your favorite Lisp implements single-floats as IEEE 754 doubles, and our library passes single-floats into IEEE 754 single-precision LAPACK routines, you'll lose accuracy when you copy back the results from LAPACK).

Observations:


Last modified: Sun Oct 8 13:26:00 CST 2006