LLA reorganization almost done

As many of you have noticed, the Lisp Linear Algebra library is currently undergoing a major reorganization. The new branch is incomplete, but what is there is perfectly is usable and passes all unit tests. The purpose of this post is to give a preview of changes and invite LLA users to experiment with the new version before it is "released".

Why the new version?

Also relying on LAPACK, the previous version of LLA used a column-major representation. This meant that even though matrices were represented using Common Lisp arrays (in particular, vectors) under the hood, they had to be wrapped to indicate that they were column-major. However, less then a year ago the LAPACK Team finally standardized the C interface for the library, which now provides row-major matrices. This means that you can use Common Lisp arrays in LLA:

(let ((a #2A((1 2)
             (3 4)))
      (b #(5 7)))
  (solve a b))

which computes \(A^{-1}b\), evaluating to

#(-3.0d0 4.0d0)

Note that the result contains double floats instead of integers: in fact, in implementations that support it, the element type of the array will be double-float.

This is because most LLA operations work like this:

1. Determine the narrowest common element type of each argument. Unless the element type of the array is something specific that LLA can recognize (eg single and double floats and their complex versions), LLA needs to traverse the array, which costs a bit in terms of speed. Use arrays with specialized element types to avoid this cost. Integers and rationals are converted to floats (you can control which, double is the default).
2. Determine the narrowest common element type of the pertinent arguments. For example, in the code above, LLA decided to use double floats, which also determines the element type.

The use of Common Lisp arrays made LLA much more transparent without any noticable overhead, especially if you use arrays that have the appropriate element types (instead of T). Even if you don't, LLA functions will return arrays that have specialized element types, so if you chain LLA operations, your code should be efficient without any particular effort.

Hermitian and triangular matrices are now wrappers around simple two-dimensional arrays. You should use as-array to convert to a dense matrix.

Required libraries

You need an implementation of LAPACK that provides the C interface. Currently, the two available choices are LAPACKE and Intel's MKL. I have been using the latter to develop the new version of LLA, so I don't have instructions on how to compile and set up LAPACKE with other libraries. If you figure it out, please write it up and I will post it in the README file of LLA.

If you install MKL on Linux, all you have to do is tell the linker where the libraries are: for example, all I had to do was create the file /etc/ld.so.conf.d/mkl.conf with contents

/opt/intel/mkl/lib/intel64
/opt/intel/composerxe/lib/intel64

and run ldconfig.

Portability

LLA uses CFFI, and should be perfectly portable. However, implementations that don't support arrays with element types single-float, double-float, (complex single-float) and (complex double-float) will suffer a performance penalty because LLA has to check the array elements to determine the narrowest type.

State of the new version

Everything works except eigenvalues and singular value decompositions. Drop me a line if you are using the new version and need either one of those, I am happy to implement them — it should be relatively easy, but I have kept features to the minimum while the library is in flux.

In principle, LLA should be clever enough to avoid copying of arrays when possible (ie no type conversion is necessary and the implementation supports foreign arrays — most of them should), but currently I am using the CFFI backend for testing.

As always, feedback, bug reports and suggestions are very welcome.

Validate XHTML 1.0