Re: GSL 2.0 roadmap (one man's view)

["Robert G. Brown" <rgb@phy.duke.edu>]

On Tue, 15 Sep 2009, Gerard Jungman wrote:

> On Mon, 2009-09-07 at 16:06 +0100, Brian Gough wrote:
>> At Thu, 27 Aug 2009 17:15:39 -0600,
>> Gerard Jungman wrote:
>>>   The important notion of slicing is (partially) implemented in GSL
>>>   in terms of the "view" concept. One can construct submatrices as
>>>   views of given matrices, change the stride of vector data by
>>>   creating vector views, etc. But there are clear flaws in the
>>>   design. The design does not express the obvious idea that
>>>   a "view" is itself a "thing", simply because the view classes
>>>   do not have an inheritance relationship to the main classes.
>>
>> I agree that the scheme is not as elegant as it could be in other
>> languages.  The view types are forced by the nature of const in C --
>> it's not possible to place the views and vectors/matrices on an equal
>> footing and preserve constness, unfortunately.  If there is a way
>> round that, I'm not aware of it.
>
>> Originally a view was essentially a vector with another name, but
>> there was no barrier to writing expressions which discarded constness
>> without any complaint from the compiler.  To prevent that in C, the
>> type had to be "wrapped" in a struct which is why one has to write
>> &view.vector or &view.matrix to access it.
>
>
> I really don't understand this. In my head, I can see a solution
> which has nothing to do with const-ness issues. I think it would
> just work. I could type it in and we could look at it. It definitely
> involves changing the way vector and matrix are done, but I don't
> think it would be a big hairy deal.
>
>
> On Mon, 2009-09-07 at 14:21 -0400, Robert G. Brown wrote:
>>
>> The biggest issue would/will probably be rationalizing the views of
>> vector and matrix so they are sufficiently portable and easy to e.g.
>> pass in and out of ODE solvers and everything else consistently.
>
> This is right to the point. It's exactly what I'm getting at.
> There is just something brain-damaged about the way it is done now.

Bless you.  And remember, it would be really nice if there was a certain
scalability in dimension.  The world doesn't stop at *a, **b -- there is
***c, ****d, etc where a number of sciences need at least through six
dimensions and back when we were kicking this around on list three or
four years ago somebody (doing relativity?) said that they use eight or
nine.  Tensors, especially tensors that don't "have" to be square, are
very useful.

In my own code, I have used at least ****d tensors that were allocated
as a single data block and the indices repacked into the tensor and then
cast to a type.  You can bass the data block to an ODE solver as a
vector, but you can write the actual deriv routine using the tensor
pointer, which means that it is actually intelligible and debuggable
without four layers of index arithmetic wrapped up on some sort of macro
or horrorshow code.

I don't remember if I mentioned it before, but multidimensional
integration (on spaces that may or may not be flat of dimension greater
than or equal to two) would be a lovely addition to the next release as
well.  There's C++ code out there (Hint) and I did a port of one of the
Hint routines at one time, but e.g. integration over a sphere of
functions with Y_lm's in them is another extremely common problem in
math and physics, and spherical domains in spherical coordinates do not
map well into nested one-D quadrature routines...

    rgb

>
> --
> G. Jungman
>
>

Robert G. Brown	                       http://www.phy.duke.edu/~rgb/
Duke University Dept. of Physics, Box 90305
Durham, N.C. 27708-0305
Phone: 1-919-660-2567  Fax: 919-660-2525     email:rgb@phy.duke.edu