Hi Jeff,
(1) Create a wrapper for the standard BLAS library. Since an optimized BLAS is available for various architectures with the same interface, this would provide a way of optimizing the library for a particular architecture.
I started working on this a few months back, but it's actually quite a large endeavor. Well, actually, BLAS isn't too bad. It's LAPACK that's a monster. IIRC, I finished porting BLAS (http://www.isi.edu/~hdaume/HBlas)...look at "documentation" to see the layout of the libraries.
(2) Encapsulate this (using a monad?) in a way that provides safe access to the imperative routines from Haskell.
On my list of things to do, but again, haven't touched it in a while.
(3) Using this as a foundation, define a matrix algebra in a "smart" way which avoids redundant evaluations. For instance, if A is an arbitrary invertible matrix, we know that A - A = 0 A + 0 = A A * Indentity = A A * inverse A = Identity so if we write an expression like "X = A * inverse A + B - B" we should expect our library to compute the correct result, "Identity" without actually performing any inversion, multiplication or elementwise matrix addition.
This is a bit more complicated, since in general, knowing if A=A is just as expensive as doing the addition/etc. One way off the top of my head to model it on top of the current HBlas library is to associate with each matrix a unique ID and a boolean. Then, when matrices are created from scratch you generate a new ID and set the boolean to false. If you invert the matrix, you leave the ID the same, but flip the boolean. The problem is with something like "A + B - B" versus "B - B + A". One of these must actually do the addition, depending on how the associativity is defined. Otherwise, in the first case, if you do the addition inner-most, you won't know "soon enough" that you don't really need to do it...
(4) Create a Haskell module which exposes a pure functional interface, hiding the imperative subsystem from the user.
I don't think this will work. I think you need the imperative interface, otherwise you are going to have to resort to copying of the arrays to maintain purity. IMO this is the wrong thing to do.
My questions are: (a) Is this idea sound? (b) Has someone already done this? (c) Is there a better way of doing efficient linear algebra in Haskell? I would greatly appreciate any advice in this regard.
I hope I answered some questions.