Re: [Caml-list] Matrix libraries

[Markus Mottl <markus@oefai.at>]

On Sat, 13 Dec 2003, Oleg Trott wrote:
> 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)

This is true for many functions in LAPACK. E.g. there is none for
computing the inverse of a matrix directly, because this is inefficient
and not very stable if you just want to solve linear equations. You'd
have to take the detour over e.g. xgetrf+xgetri.

> 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)

Maybe this is just a misunderstanding: some functions do return matrices
(and vectors), but it's always the argument that is overwritten anyway.

There is probably no purpose to having these functions return anything
if a non-optional argument containing the result is overwritten. This
will most likely just confuse people. I have changed this now. Thanks
for pointing this out!

In any case, LACAML is supposed to stay a low-level interface to
BLAS/LAPACK. Others are working on higher-level ones, e.g.:

  http://www.math.ucsb.edu/~lyons/camlFloat/index.html

I am probably going to meet up with one of the authors, Bill Lyons,
soon, who happens to come to a nearby conference, and we are intending
to talk about coordinating our work. It might be interesting to factor
out the low-level stuff using LACAML and leave higher-level functions
and operators to CamlFloat.

> 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)

I'd also like to see something like that, but I am skeptical that this
is realistically possible. Writing the C-code for the interface is very
little work. It's the OCaml-code around it that costs time to write,
especially error handling.

> 3.) Regarding "WORK" arguments. Why not have a shared workspace:
[snip]
> 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.

Yes, that's exactly the problem: what about SMP-machines and threading
then? It's really not this difficult to wrap all LACAML-functions that
use workspaces into ones that share a pool of them if threading is not
an issue.

> 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.

I had indeed thought about this, but that would have made it more
inconvenient to people who want to keep accessing "a" directly using
the Bigarray-module and the .{}-notation. CamlFloat, the higher-level
library, has such a matrix type (introduced by the constructor "Dense").
I suppose the developers are intending to add further ones for triangular,
band, symmetric matrices, etc.

> 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. 

This would collide with 4), wouldn't it? Vector arguments in LAPACK
don't have anything like a "leading dimension" so this wouldn't help you
unless you allow copying data. If, however, we keep using bigarrays as
concrete representation for matrices, then you could as well just call
"genarray_of_array2" followed by "reshape_1" in the Bigarray-module.

This is a bit inconvenient, I agree, but feel free to complain to INRIA
to introduce some more direct conversions in the Bigarray-module :-)

> 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 :-)

Ok, I have relaxed the conditions of the license. Please contribute
now :-)

Regards,
Markus

-- 
Markus Mottl          http://www.oefai.at/~markus          markus@oefai.at

-------------------
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