[racket] math/matrix

Neil Toronto writes:

 > One thing we should really do is get your LAPACK FFI into the math 
 > library and have `flmatrix-solve` use that, but fail over to Racket code 
 > systems that don't have LAPACK. If I remember right, it would have to 
 > transpose the data because LAPACK is column-major.

Sounds good. Is that LAPACK FFI already available somewhere?

 > Some thoughts, in no particular order:
 > 
 >   1. Because of transposition and FFI overhead, there's a matrix size 
 > threshold under which we ideally should use the code below, even on 
 > systems with LAPACK installed.

Transposition can be avoided in many practically relevant situations
(e.g. symmetric matrices), so such decisions should be taken on a
case-by-case basis.

As for the FFI overhead, even after reading the introduction to the FFI
documentation I have no clear idea of how important it is.  It seems
(but I am not at all certain) that there are different vector-like
data structures, some of which are optimized for access from Racket and
others for access via C pointers. If that's true, it may be of interest
to have a constructor for arrays and matrices that live in "C-optimized"
space. Or even in "Fortran-optimized" space, using column-major storage.

 >   2. Because of small differences in how it finds pivots, LAPACK's 
 > solver can return slightly different results. Should we worry about that 
 > at all?

I'd say documenting the issue is sufficient.

 >   3. A design decision: if a matrix contains just one flonum, should we 
 > convert it to (Matrix Flonum) and solve it quickly with 
 > `flmatrix-solve`, or use the current `matrix-solve` to preserve some of 
 > its exactness?
 > 
 > I lean toward regarding a matrix with one flonum as a flonum matrix. 
 > It's definitely easier to write library code for, and would make it 
 > easier for users to predict when a result entry will be exact. 

I agree. I'd like to see a clear use case for any other approach.
Predictability is  important.

Konrad.