Re: [Caml-list] Matrix libraries

[Oleg Trott <oleg_trott@columbia.edu>]

On Saturday 13 December 2003 11:36 am, Markus Mottl wrote:
> On Sat, 13 Dec 2003, romildo@uber.com.br wrote:
> > I am looking for a good library for numerical
> > matrixes manipulation in OCaml. Can anybody help me?
>
> Although it has already been on the web for a few weeks, I hadn't actually
> announced it yet, waiting for comments of some co-developers:
>
> Lacaml is available in a new major version. This library interfaces
> the Fortran libraries BLAS and LAPACK for heavy-weight linear algebra
> (i.e. matrix computations). New features include support for complex
> transformations (complex numbers), and a convenient way of accessing
> submatrices using labels. As usual, you can choose either single or
> double precision computations. The computations can run in parallel
> on SMP-machines.
>
> You can download the library here:
>
>   http://www.oefai.at/~markus/home/ocaml_sources.html#LACAML
>
> I'd be grateful for feedback!

All right! I was actually thinking of giving it lately. First of all, I think 
the library is very interesting. Here are some of my assorted thoughts 
regarding improving it.

1.) I think a matrix library has to have >= 2 levels of user-frienliness. One 
level should be easy to use and have a MATLAB feel to it, i.e. if someone
wants the eigenvalues of a real matrix, he should just be able to write
    e_vals = eigenvalues(A)
or
   (e_vals, e_vecs) = eigensystem(A)
without having to worry about details. This should be built on top of a 
low-level interface similar to LAPACK functions themselves. E.g. DGEEV does 
not even return the eigenvectors corresponding to complex eigenvalues, but 
rather something from which they can be reconstructed (which is what
one might actually want sometimes)

Right now, it appears to me that LaCaml is close to the low level, but not 
quite there (because e.g. some LAPACK interfaces _return_ a matrix rather
then letting the user modify some piece of memory in place like in FORTRAN)

2.) Currently, LaCaml interfaces only a fraction of LAPACK. LAPACK is big, and 
interfacing all of it by hand is a huge job. On the other hand, MOST of the 
information needed for the low-level interface can be parsed from function 
declarations in *.f files (I kept thinking about this while using LAPACK in 
my C++ programs). Maybe the way to interface all of LAPACK is
to parse those *.f files into some intermediate form, let humans edit it (e.g. 
add minimum workspace size expressions, etc. - obvious from docs, but hard to 
parse/NLP), and then auto-generate the whole interface from the result. This 
will have the addtional benefit of being consistent (If you just let a bunch 
of developers contribute interface code directly, they will bring their own 
"style" to their parts of the interface)

3.) Regarding "WORK" arguments. Why not have a shared workspace:

(* module Workspace, untested *)

(* private *)
let workspace = ref (Array1.create float64 100)

(* public *)
let size () = Array1.dim !workspace

let resize n = workspace := Array1.create float64 n

let delete () = resize 0

let get_float64 n = if size () >= n then !workspace 
                                                else (resize n; !workspace)

This way, one does not need to "hoist" workspace allocation out of loops 
manually. OTOH performance does not suffer from allocating workspace
for each function call. I don't know how this would inter-operate with 
threading or SMP though.

4) Submatrices. For a function f that takes a matrix argument "a" of size m by 
n, you normally use something like

f ar ac a m n

where ar and ac refer to the "beginning" of the submatrix. Even though they 
can default to 1, this is kind of inconvenient compared to

f a

I would propose using the latter everywhere, and letting the "mat" type 
include ar, ac, m and n.

submatrix : int -> int -> int -> int -> mat -> mat
(* submatrix x1 y1 x2 y2 a *)

would provide "views" into matrix a.

5) I think a function that lets one view matrices as vectors vec_of_mat is 
needed. They are all just ordered sets of numbers after all. 

6) License. Have you thought of adding a static linking exception to your LGPL 
license? Even INRIA did with Caml standard library. I might even consider 
contributing then :-)

-- 
Oleg Trott <oleg_trott@columbia.edu>

-------------------
To unsubscribe, mail caml-list-request@inria.fr Archives: http://caml.inria.fr
Bug reports: http://caml.inria.fr/bin/caml-bugs FAQ: http://caml.inria.fr/FAQ/
Beginner's list: http://groups.yahoo.com/group/ocaml_beginners