This manual is for Lisplab version 0.1, updated 10. September 2009.

Copyright © 2009 Joern Inge Vestgaarden

Lisplab manual
**************

Introduction
************

Lisplab is a mathematics library in Common Lisp. It is placed under the Gnu general public license and offers an easy-to-use and rich programming framework for mathematics, including linear algebra, Fast Fourier Transform, special functions, Runge-Kutta solver, infix notation, and general matrix utility functions. The name Lisplab is inspired by Matlab and Lisplab offers much of the same kind of programming style as Matlab, with high level manipulation of matrices. Contrary to Matlab, Lisplab can do a lot more than just matrix manipulation because Lisplab is based on Common Lisp. Hence, you get all the benefits of lexical scope, dynamic scope, macros, CLOS, fast execution, and working in an free general purpose programming language. And best of all: you can enjoy you favorite data-types in addition to the matrices: functions, hash tables, structures, classes, arbitrary precision  integers, rationals, and lists.

Lisplab is not unique in building Matlab-like syntax on top of Common Lisp. Other Common Lisp matrix libraries are Matlisp, Femlisp and NLISP. Lisplab itself was started as a branch of Matlisp, but has now moved far from the original code mass.

Getting started
***************

Dependencies
============

In order to enjoy the full power of Lisplab you must install some foreign libraries. These are

  * BLAS -- Basic Linear Algebra Subprograms. Preferably the Atlas build.

  * LAPACK -- Matrix library. Preferably the Atlas build.

  * FFTW -- The fastest Fast Fourier Transform available.

Lisplab has been developed with SBCL, SLIME and ASDF on Linux, and there are yet unnecessary bindings to these platforms.

  * Some of the optimized lisp code uses the SBCL macro truly-the. This must be dealt with.

  * The FFTW FFI is only for SBCL.

  * The Matlisp FFI should be portable to other lisps and Windows, but it has not been tested.

  * The *READ-DEFAULT-FLOAT-FORMAT* must be double-float.

Except from this, Lisplab should be self-contained and not depend on any other projects.

Installing
==========

Lisplab is ASDF installable, but before you come so far you need to specify the location of the foreign libraries. You specify these in three special variables,

  * *lisplab-libblas-path*

  * *lisplab-liblapack-path*

  * *lisplab-libfftw-path*

that live their lives in the Common-Lisp-User package. You can either assign them on the top-level, in you Common Lisp installation file, in  start.lisp, or probably the easiest: directly in lisplab.asd.

ASDF sub-systems:

  * Lisplab -- the full Lisplab installation.

  * Lisplab-base -- the part of Lisplab without external dependencies.

  * Lisplab-matlisp -- FFIs to BLAS and LAPACK. These are modified version from Matlisp.

  * Lisplab-fftw -- FFI to FFTW for Fast Fourier Transform.

  * Slatec -- special functions, generated from Fortran by f2cl. Originally made for Maxima.

  * Quadpack -- integration routines, generated from Fortran by f2cl.

The simplest way to test Lisplab is to load start.lisp. If you have problems loading, first look at start.lisp and see if you can hack it. Then look at lisplab.asd.

To install Blas, Lapack, and FFTW, if you are lazy, and lucky enough to administer a Debian or Ubuntu machine, you typically do

     # aptitude install libatlas3gf-base
     # aptitude install libfftw3-3

Naming conventions
==================

  * The matrix classes and constructors follow the naming convention from BLAS where you give names based on element type and matrix structure. The most used types are f - float, d - double, c - complex float, z - complex double float, while for matrix structure ge - general,  di - diagonal, and many more. So matrix-dge is a general matrix with double float elements, while matrix-zge is a general matrix with complex double float elements.

  * The generic functions of the basic algebra start with a dot: .+, .-, .*, ./, .^, .sin .cos, .tan, .besj, .re, etc. On numbers these functions work as the non-dotted Common Lisp functions and on matrices they work element-vise.

  * Linear algebra functions tend to start with m: m*, minv, mmax, mtp, etc., but this conventions is not strictly followed.

  * The naming convention of files follow the layered structure of Lisplab, with level0 to level3.

  * See Structure, for more about the naming conventions of matrix classes.

Status - past and future
========================

Lisplab contains a lot of linear algebra stuff, but in the future it is hope it can be broad mathematics programming environment, not just linear algebra.

Lisplab has been developed for physics simulations and data handling. Lisplab started as a refactoring of Matlisp, but the code has been to large degree rewritten, except the interfaces to Blas and Lapack. Currently, Lisplab and Matlisp have more or less the same functionality. Lisplab differ from Matlisp in the following ways

  * Implementation.

  * Layered structure (see Structure.)

  * Shorter names.

  * Lisplab uses the standard Lapack and Blas libraries, not a special build.

  * Rich matrix class hierarchy.

The future I will mainly do minor changes and bug-fixes, since it now covers my basic needs. I will only add new modules when I personally needed them.

However, there are many exiting extensions that can be made, such as

  * Parallel computation, e.g. using MPI.

  * More native linear algebra routines, e.g. eigenvalue computation.

  * New non-matrix algebra items, e.g. quaternions, polynomes or arbitrary precision floats.

  * Symbolic manipulation. Could make something like ginac in c++.

  * New matrix optimization for new usage, e.g. integer matrices for image processing or cryptography.

  * Interface to new foreign libraries, e.g. GSL.

So it this sounds interesting, please contact if you want to contribute.

Bugs and limitations
====================

The purpose of Lisplab is to be a platform for mathematical computations. From this perspective it is clear it will never be complete. Also, since there is no spec it is not obvious what is a bug and what is not!

Hence, the list in this section must be read as non-systematic gathering of problem features.

  * Lisplab runs only on SBCL. Lisplab is mainly ANSI Common Lisp, so just minor changes in the build-system should make it run on other lisps, but the problem is that it will most probably be slow. It should be fast on CMUCL, though.

  * Lacks a formal spec.

  * Poorly tested.

  * Lacks error checks (but these should not be made before a spec!)

Missing features:

  * There is no way to iterate through the elements of a general matrix in a fast way. (The map functions are currently the only thing, but these are structure agnostic and also not fast.) There should maybe be an macro w/matrix.

  * There should be linear algebra primitives, like row exchange, in level 2, so that level 3 can be made entirely without knowledge about internal structure of matrices. (Structure similar to blas - lapack)

  * Integer matrices.

  * Vectorized execution of operations.

  * Numerical integration.

  * Symbolic math. Should be separate module, only with knowledge of the dotted algebra generic functions.

  * The dotted algebra should also work on functions so the one could write (.+ (lambda (x) (+ x 1)) 3) and get a new functions as result. It might even be possible to make beautiful optimizations this way.

Tutorial
********

Starting
========

Currently Lisplab only runs on SBCL, and there is a small porting job (mainly the FFIs and a few uses of SB-EXT package) to get it run on other platforms. If you have SBCL you must make sure that lisplab.asd is (or pointed to) in asdf:*centeral-registry*. Then type

     CL-USER> (require :lisplab)
     CL-USER> (in-package :ll)

The main package of Lisplab is lisplab, with nickname ll. The tutorial assues that you are in the  lisplab package or use it. There is also a package lisplab-user, with nickname ll-user that uses lisplab.

Matrix classes
==============

The matrix classes hierarchy has three lines of inheritance. The first is on structure, the second is on element type, and the third is on implementation. This hierarchy is similar to the stream example in Object-Oriented Programming in Common Lisp, by Sonya E. Keene.

The most common classes are matrix-dge which is a general, double-float matrix with any implementation and matrix-zge which is a general, complex double-float matrix with any implementation.

The best way to inspect the class structure is to play with the low level constructor make-matrix-instance, for example,

     LL> (make-matrix-instance '(:z :ge :any) '(2 2) 1)
     #<MATRIX-ZGE  2x2
     #C(1.0 0.0) #C(1.0 0.0)
     #C(1.0 0.0) #C(1.0 0.0)
      {B529741}>

where types are :d, :z, or :any, structure is :ge or :d, and implementation is :any, :lisp, :blas. Not all combinations have a class, since not all combinations are usefull. More types and structures will be added in the future.

Note that with these three lines of inheritance there are potentially incredibly many potential class. Let's say we have 10 structures, 10 types and 4 implementations, giving 400 classes. Not all of these are useful and only few are created by default. To add new classes in a structured way, see Structure.

Matrix construction
===================

The constructor make-matrix-instance is not the primary choice to create a matrix. Its better to use a special constructor. Try

     LL> (dnew 0 2 2)
     #<MATRIX-DGE  2x2
     0.0 0.0
     0.0 0.0
      {AD78C91}>

     LL>(drow 1 2)
     #<MATRIX-DGE  1x2
     1.0 2.0
      {AD99AD9}>

     LL> (dcol 1 2)
     #<MATRIX-DGE  2x1
     1.0
     2.0
      {AEF4299}>

     LL> (dmat (1 2) (3 4))
     #<MATRIX-DGE  2x2
     1.0 2.0
     3.0 4.0
      {AF1F9A9}>

where dnew, dcol, and drow are functions, while dmat is a macro. Similarly, there are znew, zcol, zrow, and zmat for double float matrices and mnew, mcol, mrow, and mmat for any matrices. The latter take matrix class as first argument.

Often you want to create a matrix of the same type as a input matrix. Then you can use mcreate. It's mainly useful when creating methods that should operate on many matrix types.

To create matrices from something else, use convert

     LL> (convert '((1 2) (3 4)) '(:d :ge :any))
     #<MATRIX-DGE  2x2
     1.0 2.0
     3.0 4.0
      {B0998B1}>

Convert also converts between matrix types.  If the matrix contents cannot be converted directly (e.g., conversion from complex to real), use copy-contents instead.

Matrix element reference (mref and vref)
========================================

To access a matrix element use the generic function mref

     LL> (mref (dmat (1 2) (3 4)) 0 1)
     2.0

Matrix references are zero based. The mref is settable.

Furthermore, all matrices can also list their elements as vectors (similar to row-major-aref for arrays). To access a matrix as vector use the generic function vref

     LL> (vref (dmat (1 2) (3 4)) 1)
     3.0

which is also settable. Note that the matrices are column major order.

The dotted algebra (.+, .*, .sin, ...)
======================================

Central in Lisplab is an algebra (in the widest possible sense of the word) with the functions .+, .-, .*, ./, and .^. These are generalization of +, -, *, /, and ^. For numbers they work the same,

     LL> (.+ 1 2)
     3

The functions .+, .-, .*, ./, and .^ are mainly wrappers that call reduce on its arguments with the generic functions .add, .sub, .mul, .div, and .expt. If you want to extend the algebra, e.g. with polynomials, you should add specializing methods to these generic functions.

The element-vise functions by conventions have names starting with a dot (period). The dotted algebra works on the single elements and is structure agnostic. The functions of the dotted algebra introduce a programming style much the same as for Matlab, and combined with the linear algebra functions, you can write compact programs, without the slow single element references.

     LL> (.sin (drow 0 1))
     #<MATRIX-DGE  1x2
     0.0 0.8414709848078965
      {BB1FE21}>

And so also for the algebraic operations

     LL> (let ((a (drow 0 1))
               (b (dcol 0.1 0.2)))
            (.+ a (.* b 2)))
     #<MATRIX-DGE  1x2
     0.2 1.4
      {BC5AC51}>

Note here that the .+, .*, etc. ignore the structure so that it is completely valid to add a row and a column matrix. Note also that the scalar is multiplied element-vise.

The dotted algebra is to some extent also implemented for general element matrices, so that you can work with fractions and integers rather than floats

     LL> (.+ (row 'matrix-ge 1/2 3/2) 1)
     #<MATRIX-GE  1x2
     3/2 5/2
      {B5AE9B9}>

Linear algebra
==============

The linear algebra functions (unlike the dotted algebra) feel and maybe also change the matrix structure. The linear algebra functions often start with m, although this conventions is not strictly enforced. Currently there is no a wide spectrum of linear algebra function, but you find matrix multiplication, matrix inversion, transpose, conjugate transpose and eigenvalues.

Matrix multiplication

     LL>(let ((a (dmat (1 0) (0 -1)))
              (b (dcol 0.1 0.2)))
           (m* a b))
     #<MATRIX-DGE  2x1
     0.1
     -0.2
      {BD45221}>

Matrix inversion

     LL> (minv (dmat (1 2) (-2 1)))
     #<MATRIX-DGE  2x2
     0.2 -0.4
     0.4 0.2
      {BD85161}>

Transpose

     LL> (mtp (dmat (1 2) (-2 1)))
     #<MATRIX-DGE  2x2
     1.0 -2.0
     2.0 1.0
      {BD99B11}>

Conjugate transpose

     LL> (mct (zmat (1 (* 2 %i)) (-2 1)))
     #<MATRIX-ZGE  2x2
     #C(1.0 -0.0) #C(-2.0 -0.0)
     #C(0.0 -2.0) #C(1.0 -0.0)
      {BDC1DC9}>

Eigenvalues

     LL> (eigenvalues (dmat (1 2) (0 0.5)))
     #<MATRIX-DGE  2x1
     1.0
     0.5
      {B2C5B19}>

Eigenvectors

     LL> (eigenvectors (dmat (1 2) (0 0.5)))
     (#<MATRIX-DGE  2x1
     1.0
     0.5
      {B3368C1}>
      #<MATRIX-DGE  2x2
     1.0 -0.9701425001453319
     0.0 0.24253562503633297
      {B3361A9}>)

Some of the linear algebra functions also work for general element matrices

     LL> (let ((a (mat 'matrix-ge (1 0) (0 -1)))
               (b (col 'matrix-ge  1/2 2/3)))
           (m* a b))
     #<MATRIX-GE  2x1
     1/2
     -2/3
      {B5DCC79}>

     LL> (minv (mat 'matrix-ge (1 2)(-2 1)))
     #<MATRIX-GE  2x2
     1/5 -2/5
     2/5 1/5
      {BF68821}>

Matrix IO
=========

Delimited files are read with dlmread. Delimited files are written with dlmwrite. You can also write matrices as images with pgmwrite and pswrite, writing portable graymap and postscript respectively.

Note that the IO functions are currently in a poor state.

Matrices without store
======================

The class function-matrix implements matrices with functions and has no store.

The contents can for example be given directly by a rule, using the macro funmat

     LL> (funmat '(2 2) (i j) (if (= i j) 1 0))
     #<FUNCTION-MATRIX  2x2
     1 0
     0 1
      {AC0F489}>

The funmat (and its generalization fmat that creates any matrix type) it also useful to create grids.

Function matrices are also used to view a part or restructured other matrix with view-matrix, view-col, or view-row.

Matrix map
==========

There are two methods for mapping of matrices: mmap and mmap-into. For example

     LL> (mmap 'matrix-dge #'realpart (zmat (1 %i) (-%i 2)))
     #<MATRIX-DGE  2x2
     1.0 0.0
     0.0 2.0
      {A9BD489}>

The output matrix takes the structure from the first arguments, but ignores in general matrix structure. If first argument t output type is same as type of first matrix argument.

Ordinary functions
==================

These are: .sin, .cos, .sin, .tan, .sinh, .cosh, .sinh, .tanh, .asin, .acos, .asin, .atan, .asinh, .acosh, .asinh, .atanh, .conj, .re, .im, .abs, .sqrt, .exp.

Special functions
=================

These are: .besj, .besy, .besi, .besk, .besh1, .besh2, .ai, .gamma.

Infix notation
==============

Infix input is with the macro w/infix, where you must have spaces before and after the operators

     LL> (w/infix
            (let ((x 3))
               (1 .+ 2 .* x)))
     7

The w/infix messes as little as possible with the Lisp semantics, so that if you have a lot of formulas just wrap all of it inside the macro. The infix math also works with the functions +, -, *, / and ^.

Structure
*********

Design principles
=================

Design principles for the full library

  * Lisplab is free software.

  * It makes a homogeneous platform for all kinds of mathematical calculations. (So it's a lot more than just a matrix library)

  * User applications should need to stay only in Common Lisp. (There should be no need for optimized math in FFIs or special languages like Maxima)

  * Every common mathematical operator and function is represented by a CLOS generic function. This is called the dotted algebra.

  * Modular structure (Inspired by GSL).

  * Trust the Lisp system and use foreign code as little as possible.

  * Avoid programming mathematical algorithms in macros. Despite the advantages (fast and generic at the same time) it is hard to understand and debug.

  * Error checks is primarily callers responsibility, not Lisplab's!

  * To steal as much code as possible from as many as possible (I love free software).

Design principles for the matrix part

  * Layered structure where dependencies are primarily on the layer below -- not vertical within the layer.

  * Layer 0 is mainly interfaces and generic functions, not implementations. Level 1 and 2 are small. This structure should encourage modularity.

Package structure
=================

So far, there is only one main package, called, you might guess it: lisplab. Except from that there are only a few special packages for generated code and FFIs: Slatec, Blas, and FFTW. For test code an applications you have the package lisplab-user.

The four levels, 0 -- 3.
========================

The Lisplab matrix and linear algebra code has a layered structure with four levels, 0 -- 3, where

  * Level 0 is matrix independent and contains the dotted algebra generic functions, the dotted algebra methods specialized for numbers, and helper functions and macros.

  * Level 1 defines matrices and defines and implements mref, vref,cols, rows, size and  make-matrix-instance.

  * Level 2 defines and core functionality related to matrices, such as the dotted algebra methods, matrix constructors (dnew, dcol, etc.) and other matrix helper functions, such as mmax, mmin, circ-shift, etc.

  * Level 3 is everything else that uses matrices, including linear algebra, FFTs, solvers, etc.

The levels are unequal in size: level 0 is fairly large, level 1 is extremely small, level 2 is fairly small, and level 3 is large and I hope it will never stop to grow.

The intention with the structure are the following

  * To minimize the work needed to add new matrices implementations. Actually, you need only to implement level 1, provided you inherit from matrix base. Otherwise you need also to implement level 2.

  * To encourage code reuse. For example, if you add a new type to the dotted algebra, such as polynomials, symbolic expressions, or arbitrary precision floats, you can immediately perform matrix operations on them without write one single extra line of code.

  * To manage optimizations and special implementations by overloading for specialized matrix types, mainly at level 2, but also at level 3.

  * To make foreign libraries available in Common Lisp.

  * To create a closed, constant and extensible system of operations.

The filenames often also denote the levels. If a filename does not denote the level, it is most probably level 3, or outside the level system (non matrix code).

Matrix class hierarchy
======================

All matrices are subclasses of matrix-base, and as far as possible the generic functions specialize on one or many of its subclasses.

The matrix class hierarchy has three independent lines of inheritance

  * On structure

  * On element type

  * On implementation

The structure is inspired by the stream example in Object-Oriented Programming in Common Lisp, by Sonya E. Keene.

The structure
-------------

The structure tells what kind of symmetries or other special properties the matrix has. Two matrix classes

  * matrix-structure-general, where the matrix has no known symmetries.

  * matrix-structure-diagonal, where only diagonal elements are non-zero.

  * Other types to come.

The element type
----------------

The element type classes has no other purpose than to be represents a Common Lisp types.

  * matrix-element-base, represents t

  * matrix-element-double-float, represents double-float

  * matrix-element-complex-double-float, represents complex double-float

The implementation
------------------

Since Lisplab has many competing implementation of the same generic functions, the implementation class structure tells which one too choose. There are currently four classes in a straight line of inheritance

  * matrix-implementation-base

  * matrix-implementation-lisp, use native Common Lisp if possible.

  * matrix-implementation-blas, use foreign libraries if possible.

This standard method dispatch ensures foreign library methods are chosen before native lisp. It is the responsibility of the methods to call next method if the library is not loaded.

Discussion
**********

The foreign function interfaces
===============================

The elements of Lisplab typed matrices (double-float and complex double-float) are stored as 1D simple arrays, and most Lisps will then have a SAP inside the array pointer which is binary compatible with Fortran. This adds some overhead to matrix element references, but simplifies the garbage collections.

Matlisp has an extra layer with Fortran compatibility above the ordinary C FFI. This layer is mainly unchanged compared to Matlisp, but only the SBCL FFI is in the code. The FFIs for other lisps could easily be added, but not without testing.

The FFI for FFTW is a mock-up for SBCL, and it only is for standard complex transforms and inverse transforms. Of course, a lot more can be done, but this simple version works quite OK. In long term the SBCL FFI should be replaced with a general one.

Parallel execution
==================

To request usage of more threads, call init-threads. Currently only the calls to FFTW will react to this.

Symbolic calculations
=====================

There is draft code for symbolic calculations,  but it is not yet included in the build. The symbolic code need not know anything about the rest of Lisplab, except for the generic methods of the dotted algebra (The level 0).