From mboxrd@z Thu Jan 1 00:00:00 1970 Received: (from majordomo@localhost) by pauillac.inria.fr (8.7.6/8.7.3) id QAA17442; Sun, 14 Dec 2003 16:02:00 +0100 (MET) X-Authentication-Warning: pauillac.inria.fr: majordomo set sender to owner-caml-list@pauillac.inria.fr using -f Received: from nez-perce.inria.fr (nez-perce.inria.fr [192.93.2.78]) by pauillac.inria.fr (8.7.6/8.7.3) with ESMTP id QAA17488 for ; Sun, 14 Dec 2003 16:01:58 +0100 (MET) Received: from fichte.ai.univie.ac.at (fichte.ai.univie.ac.at [131.130.174.156]) by nez-perce.inria.fr (8.11.1/8.11.1) with ESMTP id hBEF1v104655 for ; Sun, 14 Dec 2003 16:01:57 +0100 (MET) Received: from fichte.ai.univie.ac.at (markus@localhost [127.0.0.1]) by fichte.ai.univie.ac.at (8.12.3/8.12.3/Debian-6.6) with ESMTP id hBEF1rHn014305; Sun, 14 Dec 2003 16:01:53 +0100 Received: (from markus@localhost) by fichte.ai.univie.ac.at (8.12.3/8.12.3/Debian-6.6) id hBEF1qAQ014302; Sun, 14 Dec 2003 16:01:52 +0100 Date: Sun, 14 Dec 2003 16:01:52 +0100 From: Markus Mottl To: Oleg Trott Cc: romildo@uber.com.br, OCaml Subject: Re: [Caml-list] Matrix libraries Message-ID: <20031214150152.GB4945@fichte.ai.univie.ac.at> Mail-Followup-To: Oleg Trott , romildo@uber.com.br, OCaml References: <20031213135142.GA8519@gentoo.malaquias.no-ip.org> <20031213163610.GA4004@fichte.ai.univie.ac.at> <200312131638.01256.oleg_trott@columbia.edu> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <200312131638.01256.oleg_trott@columbia.edu> User-Agent: Mutt/1.4.1i X-Loop: caml-list@inria.fr X-Spam: no; 0.00; caml-list:01 oleg:01 matlab:01 vals:01 vals:01 low-level:01 lapack:01 lapack:01 lacaml:01 vectors:01 overwritten:01 non-optional:01 overwritten:01 lacaml:01 low-level:01 Sender: owner-caml-list@pauillac.inria.fr Precedence: bulk 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