[Camllist] Matrix libraries
[
Home
]
[ Index:
by date

by threads
]
[ Message by date: previous  next ] [ Message in thread: previous  next ] [ Thread: previous  next ]
[ Message by date: previous  next ] [ Message in thread: previous  next ] [ Thread: previous  next ]
Date:   (:) 
From:  Oleg Trott <oleg_trott@c...> 
Subject:  Re: [Camllist] Matrix libraries 
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 codevelopers: > > Lacaml is available in a new major version. This library interfaces > the Fortran libraries BLAS and LAPACK for heavyweight 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 SMPmachines. > > 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 userfrienliness. 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 lowlevel 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 lowlevel 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 autogenerate 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 interoperate 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 camllistrequest@inria.fr Archives: http://caml.inria.fr Bug reports: http://caml.inria.fr/bin/camlbugs FAQ: http://caml.inria.fr/FAQ/ Beginner's list: http://groups.yahoo.com/group/ocaml_beginners