gsl container designs

gsl container designs

[Gerard Jungman]

[-- Attachment #1: Type: text/plain, Size: 266 bytes --]


Here are header files for a couple different approaches to containers.
I didn't bother with any implementations; it seems obvious how to
implement most of these functions.

The designs are not complete, but they express most of
the important stuff.

--
G. Jungman


[-- Attachment #2: gsl-container-designs.tar.gz --]
[-- Type: application/x-compressed-tar, Size: 11583 bytes --]

Re: gsl container designs

[Brian Gough]

At Tue, 24 Nov 2009 17:54:46 -0700,
Gerard Jungman wrote:
> Here are header files for a couple different approaches to containers.
> I didn't bother with any implementations; it seems obvious how to
> implement most of these functions.
> 
> The designs are not complete, but they express most of
> the important stuff.

Thanks, I will study them this week.

Re: gsl container designs

[Brian Gough]

At Tue, 24 Nov 2009 17:54:46 -0700,
Gerard Jungman wrote:
> Here are header files for a couple different approaches to containers.
> I didn't bother with any implementations; it seems obvious how to
> implement most of these functions.
> 
> The designs are not complete, but they express most of
> the important stuff.

Thanks for the document, I have studied the designs this week.  It
seems that changing to design 1 / 1u / 2 would be trading one set of
problems for another.  Looking at each case, the change doesn't seem
sufficient to justify the compatibility cost.

Regarding the introductory points.

2a+b) While we may have had some kind of goal of supporting C++ I
don't think it's worth encouraging that today as the gap between the
languages has widened so much since the start of the project.

3) Non-levelised types.  These seem to be the price for type safety.
In terms of the look/feel, expressions like &row.vector and
&column.vector don't seem too unnatural to me.

4) I think this can be fixed in the current framework.

5) Since row-major was in the original design it's too much of a break
to change that.

I think the most realistic way forward is to add a multiarray type,
with rank-1, rank-2 and rank-N versions, in the existing framework and
use non-levelised types for const/non-const views to preserve
type-safety.

-- 
Brian Gough

using GSL with C++ (was Re: gsl container designs)

[James Amundson]

On 12/04/2009 12:36 PM, Brian Gough wrote:
> At Tue, 24 Nov 2009 17:54:46 -0700,
> Gerard Jungman wrote:
>    
>> Here are header files for a couple different approaches to containers.
> Regarding the introductory points.
>
> 2a+b) While we may have had some kind of goal of supporting C++ I
> don't think it's worth encouraging that today as the gap between the
> languages has widened so much since the start of the project.
>    

I, for one, use GSL almost exclusively from C++. When I'm not using C++, 
I am usually using Python. How are other readers of this list using GSL? 
Are you using a pure C environment, or do you use C++ and/or other 
languages?

I have no interest in debating the relative merits of various languages; 
I am just curious to know what is happening in practice.

--Jim Amundson

Re: using GSL with C++ (was Re: gsl container designs)

[Jochen Küpper]

On 09.12.2009, at 22:02, James Amundson wrote:

> On 12/04/2009 12:36 PM, Brian Gough wrote:
>> 

>> 2a+b) While we may have had some kind of goal of supporting C++ I
>> don't think it's worth encouraging that today as the gap between the
>> languages has widened so much since the start of the project.  
> 
> I, for one, use GSL almost exclusively from C++. When I'm not using C++, I am usually using Python. How are other readers of this list using GSL? Are you using a pure C environment, or do you use C++ and/or other languages?


C++ and sometimes Python (here I more often use scipy than pygsl)

Re: using GSL with C++ (was Re: gsl container designs)

[Jari Häkkinen]

I use GSL exclusively from C++. I cannot remember when I wrote C ... 
probably 12-13 years has passed since I did my last (latest?) C 
consultancy gig.

Jari

James Amundson wrote:
> On 12/04/2009 12:36 PM, Brian Gough wrote:
>> At Tue, 24 Nov 2009 17:54:46 -0700,
>> Gerard Jungman wrote:
>>   
>>> Here are header files for a couple different approaches to containers.
>> Regarding the introductory points.
>>
>> 2a+b) While we may have had some kind of goal of supporting C++ I
>> don't think it's worth encouraging that today as the gap between the
>> languages has widened so much since the start of the project.
>>    
>
> I, for one, use GSL almost exclusively from C++. When I'm not using 
> C++, I am usually using Python. How are other readers of this list 
> using GSL? Are you using a pure C environment, or do you use C++ 
> and/or other languages?
>
> I have no interest in debating the relative merits of various 
> languages; I am just curious to know what is happening in practice.
>
> --Jim Amundson

Re: using GSL with C++ (was Re: gsl container designs)

[Brian Gough]

At Wed, 09 Dec 2009 15:02:20 -0600,
James Amundson wrote:
> I, for one, use GSL almost exclusively from C++. When I'm not using C++, 
> I am usually using Python. How are other readers of this list using GSL? 
> Are you using a pure C environment, or do you use C++ and/or other 
> languages?

What types of functions do you mainly use?  The impedance mismatch for
the "higher-order" functions, e.g. using functions arguments and void
* parameters, seems to be the problematic area to me. 

Re: using GSL with C++ (was Re: gsl container designs)

[Brian Gough]

At Thu, 10 Dec 2009 11:44:54 +0000,
Brian Gough wrote:
> What types of functions do you mainly use?  The impedance mismatch for
> the "higher-order" functions...

Specifically in C++, where there are better ways to do it.

Re: using GSL with C++ (was Re: gsl container designs)

[Robert G. Brown]

On Wed, 9 Dec 2009, James Amundson wrote:

> On 12/04/2009 12:36 PM, Brian Gough wrote:
>> At Tue, 24 Nov 2009 17:54:46 -0700,
>> Gerard Jungman wrote:
>> 
>>> Here are header files for a couple different approaches to containers.
>> Regarding the introductory points.
>> 
>> 2a+b) While we may have had some kind of goal of supporting C++ I
>> don't think it's worth encouraging that today as the gap between the
>> languages has widened so much since the start of the project.
>> 
>
> I, for one, use GSL almost exclusively from C++. When I'm not using C++, I am 
> usually using Python. How are other readers of this list using GSL? Are you 
> using a pure C environment, or do you use C++ and/or other languages?
>
> I have no interest in debating the relative merits of various languages; I am 
> just curious to know what is happening in practice.

I use C only.  Consequently I come down in favor of using casts to type
e.g. generalized tensor forms, because pointer-based indirection can
lead to reasonably efficient, highly readable code (and because the GSL
really needs a generalized tensor form -- the computational universe of
interest to many scientists is not just one or two dimensional,
rectangular, with indices starting at 0 or 1.

    rgb

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

Re: using GSL with C++ (was Re: gsl container designs)

[James Bergstra]

On Thu, Dec 10, 2009 at 8:42 AM, Robert G. Brown <rgb@phy.duke.edu> >
I use C only.  Consequently I come down in favor of using casts to
type
> e.g. generalized tensor forms, because pointer-based indirection can
> lead to reasonably efficient, highly readable code (and because the GSL
> really needs a generalized tensor form -- the computational universe of
> interest to many scientists is not just one or two dimensional,
> rectangular, with indices starting at 0 or 1.
>
>   rgb
+1

I mainly use python, numpy, & scipy now. But one of the main reasons
for switching from gsl (aside from the python language) was natural
support and syntax that numpy provides for n-dimensional arrays.  I
typically use arrays with n <= 5.

James
-- 
http://www-etud.iro.umontreal.ca/~bergstrj

Re: using GSL with C++ (was Re: gsl container designs)

[Kevin H. Hobbs]

[-- Attachment #1: Type: text/plain, Size: 382 bytes --]

On 12/09/2009 04:02 PM, James Amundson wrote:
> How are other readers of this list using GSL? Are you using a pure C
> environment, or do you use C++ and/or other languages?

I use GSL almost exclusively from C.

I'm usually running lots of ODE simulations where the parameters are
tuned by an MPI based evolutionary algorithm.

From time to time I'll use GSL from C++.


[-- Attachment #2: OpenPGP digital signature --]
[-- Type: application/pgp-signature, Size: 252 bytes --]

Re: gsl container designs

[Tuomo Keskitalo]

Hi,

On 12/04/2009 08:36 PM, Brian Gough wrote:

> At Tue, 24 Nov 2009 17:54:46 -0700,
> Gerard Jungman wrote:
>> Here are header files for a couple different approaches to containers.
>> I didn't bother with any implementations; it seems obvious how to
>> implement most of these functions.
>>
>> The designs are not complete, but they express most of
>> the important stuff.
> 
> Thanks for the document, I have studied the designs this week.  It
> seems that changing to design 1 / 1u / 2 would be trading one set of
> problems for another.  Looking at each case, the change doesn't seem
> sufficient to justify the compatibility cost.

Do you mean the compatibility to GSL 1 types by the compatibility cost? 
When talking about GSL 2 I don't think we should give too much value to 
maintaining backwards compatibility. GSL 1 is not going to cease to 
exist, and people who have tied themselves deeply to GSL 1 data 
structures can continue to use it.

> 3) Non-levelised types.  These seem to be the price for type safety.
> In terms of the look/feel, expressions like &row.vector and
> &column.vector don't seem too unnatural to me.

Here is a crazy idea: If we take Gerards design 1, would it be too 
absurd to use a horrible big wrapper struct like

typedef struct {
   gsl_marray *m;
   gsl_const_marray *cm;
   gsl_marray_1 *m1;
   gsl_const_marray_1 *cm1;
   gsl_marray_2 *m2;
   gsl_const_marray_2 *cm2;
   gsl_marray_3 *m3;
   gsl_const_marray_3 *cm3;
   gsl_vector *vec;
   gsl_const_vector *cvec;
   gsl_matrix *mat;
   gsl_const_matrix *cmat;
   ...
} gsl_container;

and make everything a gsl_container, and then always use them like
&a.mat or &b.vec? Maybe the const types could be in another 
gsl_const_container struct..

One downside is that this would effectively mask the type of a or b in 
program code (gsl_container a; // is it a vector, matrix, or marray?), 
which might make reading code a pain..

-- 
Tuomo.Keskitalo@iki.fi
http://iki.fi/tuomo.keskitalo

Re: gsl container designs

[Robert G. Brown]

On Wed, 6 Jan 2010, Tuomo Keskitalo wrote:

> Here is a crazy idea: If we take Gerards design 1, would it be too absurd to 
> use a horrible big wrapper struct like
>
> typedef struct {
>  gsl_marray *m;
>  gsl_const_marray *cm;
>  gsl_marray_1 *m1;
>  gsl_const_marray_1 *cm1;
>  gsl_marray_2 *m2;
>  gsl_const_marray_2 *cm2;
>  gsl_marray_3 *m3;
>  gsl_const_marray_3 *cm3;
>  gsl_vector *vec;
>  gsl_const_vector *cvec;
>  gsl_matrix *mat;
>  gsl_const_matrix *cmat;
>  ...
> } gsl_container;
>
> and make everything a gsl_container, and then always use them like
> &a.mat or &b.vec? Maybe the const types could be in another 
> gsl_const_container struct..
>
> One downside is that this would effectively mask the type of a or b in 
> program code (gsl_container a; // is it a vector, matrix, or marray?), which 
> might make reading code a pain..

IMO the primary reason for "fixing" vector/matrix objects in the GSL is
to enable the single vector that represents the actual data in the
object to be presented to the user in a matrix-intuitive way, even at
some cost in code efficiency (and to try to minimize that cost or
provide further benefits to counterbalance it by enabling efficient
linear algebra routines to be matched to it, for example).

Any of us who have coded for a decade or two know perfectly well how to
address a vector via explicit displacement arithmetic, burying horrible
constructs in one's code to compute the displacement of an i,j,k element
in a vector containin aa "3d matrix".  Or, one can use C pointer
constructs that are equally unreadable for the same general purpose.  We
have also had to DEBUG this sort of thing, which is extremely
unpleasant, or understand somebody else's displacement code, which is
often unpleasant to the point of being nightmarish.

What I want is extremely simple:

* A way to construct/allocate a matrix of arbitrary dimension (out to some
(un)reasonable max, e.g. 10, likely to be higher than anybody needs for
coding purposes at least at this point, and deconstruct/free the matrix
when I'm done with it, ideally with a single call in each case, ideally
for any type of variable.  I'm perfectly comfortable with requiring a
pointer and cast to establish the type, but I want to be able to use it
just like m[i][j][k]...[r][s] when all is said and done.

* The specific dimension bounds need to be user specified.  If I want to
make a 6 dimensional triple triangular array to hold e.g. 3j symbols, I
don't want to have to allocate a rectangular array most of which is
empty wasted space.  I'm not talking about sparse matrix here, just
non-rectangular arrays.  Angular momentum indices are perfect (and
common) examples.

* Access to the underlying actual data vector (non-opaque, in other
words) so it can be passed to e.g. ODE solvers and the like that can
only be sensibly written in a general way for a vector.  I want to be
able to pass the ODE solver the vector and its dimension (and perhaps
an auxiliary similar vector for the derivatives with similar packing
and indexing), but I want to be able to write the ODE deriv routine 
using the readable matrix form.

None of this is particularly difficult to accomplish.  I do it all the
time in my own code.  In the end, one can integrate sets of ODEs where
the derivatives are indexed by angular momentum l,m, where the deriviative
expression contains sums over l',m', l",m" with precomputed tables of
e.g. Gaunt numbers, and where the derivative expression is HUMAN
READABLE and directly comparable to an expression in a physics paper or
textbook.

My problem in this regard with the GSL is that so far, the matrix and
vector constructs are pretty much useless -- they are far away from the
standard way one usually does this sort of thing in C using pointers or
just allocating a vector or matrix.

One of the ADVANTAGES of C is that it permits one to work pointer magic
to improve the readability (and sometimes execution efficiency) of code.
Opaque types and set/get calls are dark evil in both regards.  I just
want to be able to do (double) I[l][m][l'][m'][l"][m"]aand/or (int)
epsilon[i][j][k] transparently and reasonably efficiently, because both
of these occur all the time in real numerical code tht I've written.  I
can do it myself and still use the GSL, and have in the past done so.
For ODE solvers this isn't much of a problem.  For linear algebra,
however, it starts being a problem.  Then issues like column-major
become important and have to be solved in a standardized way.  This
would make even my own code more efficient, as well, as most of what I do
with the I() matrix above is basically matrix multiplication or
contraction -- summing out shared indices -- and linear algebra routines
that operate on correctly allocated structures can be 2x-3x more
efficient than ones that sum across big steps in the underlying vectors.

    rgb

>
> -- 
> Tuomo.Keskitalo@iki.fi
> http://iki.fi/tuomo.keskitalo
>

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

Re: gsl container designs

[Gerard Jungman]

On Wed, 2010-01-06 at 10:47 -0500, Robert G. Brown wrote:

First, I want to say that I agree with everything you
say. These are exactly the problems that we need to solve.


> IMO the primary reason for "fixing" vector/matrix objects in the GSL is
> to enable the single vector that represents the actual data in the
> object to be presented to the user in a matrix-intuitive way, even at
> some cost in code efficiency (and to try to minimize that cost or
> provide further benefits to counterbalance it by enabling efficient
> linear algebra routines to be matched to it, for example).

In fact, I don't even see a need to sacrifice efficiency. I think 
this goal is directly achievable at essentially optimal efficiency.
My attitude about C is that it is all about lining up the bytes.
A good C design is just some syntactic layer on top of a
rudimentary layer, where-in you are careful to line up all
the bytes. After that, optimal efficiency is pretty much
automatic. C is very friendly, once you are careful to
line everything up. I hope you see my small but non-vacuous
technical meaning.


> What I want is extremely simple:
> 
> * A way to construct/allocate a matrix of arbitrary dimension (out to some
> (un)reasonable max, e.g. 10

Check.


> but I want to be able to use it
> just like m[i][j][k]...[r][s] when all is said and done.

This is the only part which I wonder about. I thought about this
before and came to the conclusion that it is not easy to get right.
Probably possible, but very tricky.

The most annoying aspect of this is that these little problems
are trivial in C++. But maybe just beyond the capability of C.


> * The specific dimension bounds need to be user specified.  If I want to
> make a 6 dimensional triple triangular array to hold e.g. 3j symbols, I
> don't want to have to allocate a rectangular array most of which is
> empty wasted space.

See the earlier comment by Tuomo Keskitalo, about banded/packed storage.
Not something I have thought about yet, but probably not hard.


> * Access to the underlying actual data vector (non-opaque, in other
> words) so it can be passed to e.g. ODE solvers and the like that can
> only be sensibly written in a general way for a vector.

Check.


> My problem in this regard with the GSL is that so far, the matrix and
> vector constructs are pretty much useless -- they are far away from the
> standard way one usually does this sort of thing in C using pointers or
> just allocating a vector or matrix.

I find them useless as well, for exactly these reasons.


> One of the ADVANTAGES of C is that it permits one to work pointer magic
> to improve the readability (and sometimes execution efficiency) of code.
> Opaque types and set/get calls are dark evil in both regards.

Right. The layer which defines the "types" should be essentially
syntactic, just there to make it easier to organize/understand
your code. The underlying implementation should be a brutally
simple pile of bytes, and it should be visible.


--
G. Jungman

Re: gsl container designs

[Robert G. Brown]

On Wed, 6 Jan 2010, Gerard Jungman wrote:

> In fact, I don't even see a need to sacrifice efficiency. I think

I recall from my benchmarking days that -- depending on compiler --
there is a small dereferencing penalty for packed matrices (vectors
packed into dereferencing **..* pointers) compared to doing the offset
arithmetic via brute force inline or via a macro.  I ran two different
implementations of stream (one with direct pointer traversal of the
vectors, one with vector offsets repacked into **matrices to do the
tests.  The size of this penalty did decrease as gcc advanced; I haven't
run the benchmark recently and don't know how large it currently is.  It
was never so large that it stopped me from using repacked pointers for
code clarity.

> this goal is directly achievable at essentially optimal efficiency.
> My attitude about C is that it is all about lining up the bytes.
> A good C design is just some syntactic layer on top of a
> rudimentary layer, where-in you are careful to line up all
> the bytes. After that, optimal efficiency is pretty much
> automatic. C is very friendly, once you are careful to
> line everything up. I hope you see my small but non-vacuous
> technical meaning.

Agreed.  And I can see no reason that modern compilers may not have
eliminated most of the dereferencing penalty, especially if one aligns
things carefully.  However, some of my beowulfish friends who work(ed)
for compiler companies argued that one of Fortran's longstanding
advantages in numerical code efficiency is simply because of the fact
that matrices in Fortran are written in stone and so compiler writers
can optimize the hell out of them.  C pointers, OTOH -- well, one can
make triangular arrays, rectangular arrays, one can start the indices at
0 or 1 or -100 or 47 as you like and run to wherever you like.  How can
anyone pre-optimize a compiler to handle all of that?

Truthfully, I can't see why they CAN'T, but I'm not a compiler writer
and empirically there has been a small penalty between C in general and
fortran and C with direct addressing and C with dereferenced addressing
as well.  IIRC there is one more layer of address resolution in indirect
dereferencing, which is the fundamental origin of the latter.

>> What I want is extremely simple:
>>
>> * A way to construct/allocate a matrix of arbitrary dimension (out to some
>> (un)reasonable max, e.g. 10
>
> Check.
>
>
>> but I want to be able to use it
>> just like m[i][j][k]...[r][s] when all is said and done.
>
> This is the only part which I wonder about. I thought about this
> before and came to the conclusion that it is not easy to get right.
> Probably possible, but very tricky.

I've written code to do so many times -- Numerical Recipes provides the
basic idea in its matrix-packing routines.  Simply allocate e.g.

double **..*m,*v;

Then malloc v long enough to hold the actual data.  Then malloc m to
hold the first index, malloc vectors long enough to hold the second
index and assign them to these indices in a first-index loop (recurse
until the next to last index) and finally pack the appropriately offset
addresses from the actual vector v into the next to last index so that
the last index actually traverses v.  In this scheme the last index is
always the linear/fast index in the actual data vector, the other
indices are all there to permit the compiler to handle the
dereferencing.

To free the matrix, one has to work backwards -- free each successive
index one malloc'd from the right back to the left (you need the
leftmost indices to know the addresses to free at each level), and
finally free m and v themselves.  All very tedious and mechanical, but
all fully automatable and debuggable as well.  The biggest problem with
NR is that each type has an individual allocation/free command, and it
only goes to 2 dimensions.  To make a single routine pair handle more
dimensions and arbitrary types is more work and requires void typing and
casting (and/or switches to handle the unique aspects of individual
tyees internally).

> The most annoying aspect of this is that these little problems
> are trivial in C++. But maybe just beyond the capability of C.

How could that be?  People write operating systems in C.  I hardly think
that allocating sensible matrices is beyond it.

I'm not trying to be argumentative or flame, by the way.  Could you be
more explicit about why packing pointers to permit the dereferencing of
an arbitrary matrix is technically difficult (as opposed to merely
tedious) in C and would be easier in C++ (which I thought discouraged
the use of pointers altogether)?

>> * The specific dimension bounds need to be user specified.  If I want to
>> make a 6 dimensional triple triangular array to hold e.g. 3j symbols, I
>> don't want to have to allocate a rectangular array most of which is
>> empty wasted space.
>
> See the earlier comment by Tuomo Keskitalo, about banded/packed storage.
> Not something I have thought about yet, but probably not hard.

I'm not sure about "banded", but packed is (as noted) simple.  It is
just a matter of picking the lengths of the successive ***, **, * length
vectors one mallocs to fill the toplevel e.g. **** pointer, plus the use
of appropriate offsets.  In particular, NR contains examples of how to
do e.g.:

   **dmatrix(l1,u1,l2,u2)

for arbitrary lower/upper indices in the two slots (as well as the
associated free) and it is straightforward to generalize the basic
algorithm to higher dimensions and allow for arbitrary data types.  I
don't think twice about allocating a vector of structs of more or less
arbitrary length (including structs that contain vectors or matrices
themselves) and dereferencing it via a packed, cast or typed, matrix.

This is one of the ADVANTAGES of C -- one can just do what one wants to
do in the most natural way via pointers, mallocs and suitable typedefs
and/or structs.  Without a safety net, admittedly, and you're
responsible for both code integrity and efficiency with little compiler
help, but OTOH you can quite easily visualize what the underlying
assembler is going to (approximately) look like, so one CAN tweak it in
reasonable ways for efficiency without resorting to actual assembler.

>> * Access to the underlying actual data vector (non-opaque, in other
>> words) so it can be passed to e.g. ODE solvers and the like that can
>> only be sensibly written in a general way for a vector.
>
> Check.
>
>
>> My problem in this regard with the GSL is that so far, the matrix and
>> vector constructs are pretty much useless -- they are far away from the
>> standard way one usually does this sort of thing in C using pointers or
>> just allocating a vector or matrix.
>
> I find them useless as well, for exactly these reasons.
>
>
>> One of the ADVANTAGES of C is that it permits one to work pointer magic
>> to improve the readability (and sometimes execution efficiency) of code.
>> Opaque types and set/get calls are dark evil in both regards.
>
> Right. The layer which defines the "types" should be essentially
> syntactic, just there to make it easier to organize/understand
> your code. The underlying implementation should be a brutally
> simple pile of bytes, and it should be visible.

Amen, my brother.  C is a THIN veneer of upper-level language
sensibilities on top of raw assembler, and its ability to work "close to
the metal" is its power.  This power should be preserved as much as
possible in the GSL, in the process of providing a human readable
interface.  I actually wrote all of this, in a more or less GSL--like
form, several years ago.  Some of the code was ugly, but it worked
charmlike, at the cost of requiring a cast of a void type to tell the
compiler how to dereference the outermost (actual) vector index.
Internally, all the dereferencing pointer arithmetic is more or less the
same for a void (universal) type, but the pointer into the actual data vector 
has to be advanced with the correct stride, and the compiler has to know 
what that stride should be.

I can probably dig this code out and post it here (or post a link to it)
if you want to look at it.  I think I went up to 9 or 10 dimensions (way
more than I needed) and did a lot of things via brute force to make the
purely mechanical recursion obvious; I planned an eventual rewrite into
a more compact form that I never got around to.

    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

Re: gsl container designs

[Rhys Ulerich]

> I recall from my benchmarking days that -- depending on compiler --
> there is a small dereferencing penalty for packed matrices (vectors
> packed into dereferencing **..* pointers) compared to doing the offset
> arithmetic via brute force inline or via a macro.
> ......
> I haven't
> run the benchmark recently and don't know how large it currently is.  It
> was never so large that it stopped me from using repacked pointers for
> code clarity..

Mostly unscientific, but worth tossing into the mix:

Using Intel 10.1 compilers on a fairly recent AMD chip, 100,000 iterations
of doing the nested pointers approach is neck-and-neck with index arithmetic
on a 10x10 double matrix.  For the 100x100 case it takes 1.3 times longer
to iterate using the nested pointers.  Work in the inner loop "compute
kernel" is
*= against a constant scalar.  Optimization flags on -O3.  I've seen similar
behavior on recent GNU compilers.

I'm happy to provide the test code if anyone's interested.

- Rhys

Re: gsl container designs

[Robert G. Brown]

[-- Attachment #1: Type: TEXT/PLAIN, Size: 1760 bytes --]

On Wed, 6 Jan 2010, Rhys Ulerich wrote:

>> I recall from my benchmarking days that -- depending on compiler --
>> there is a small dereferencing penalty for packed matrices (vectors
>> packed into dereferencing **..* pointers) compared to doing the offset
>> arithmetic via brute force inline or via a macro.
>> ......
>> I haven't
>> run the benchmark recently and don't know how large it currently is.  It
>> was never so large that it stopped me from using repacked pointers for
>> code clarity..
>
> Mostly unscientific, but worth tossing into the mix:
>
> Using Intel 10.1 compilers on a fairly recent AMD chip, 100,000 iterations
> of doing the nested pointers approach is neck-and-neck with index arithmetic
> on a 10x10 double matrix.  For the 100x100 case it takes 1.3 times longer
> to iterate using the nested pointers.  Work in the inner loop "compute
> kernel" is
> *= against a constant scalar.  Optimization flags on -O3.  I've seen similar
> behavior on recent GNU compilers.

That sounds partly like a cache effect -- 10x10 almost certainly stays
in L1, 100x100 won't fit.  My own experience is similar, although I
don't recall the multiplier being as large as 1.3 (but then, I was doing
stream and stream-like tests with very large vectors, which means that
one spends more time in a vector streaming mode and minimizes
cache-thrashing when turning corners).  And my memory could be faulty --
I'm an old guy, after all, early Alzheimers...;-)

    rgb

>
> I'm happy to provide the test code if anyone's interested.
>
> - Rhys
>

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

Re: gsl container designs

[Gerard Jungman]

On Wed, 2010-01-06 at 22:29 -0500, Robert G. Brown wrote:

> However, some of my beowulfish friends who work(ed)
> for compiler companies argued that one of Fortran's longstanding
> advantages in numerical code efficiency is simply because of the fact
> that matrices in Fortran are written in stone and so compiler writers
> can optimize the hell out of them.  C pointers, OTOH -- well,

I think the main difference is that fortran can assume no aliasing,
which allows some extra optimization that C cannot do. The new 
'restrict' keyword is supposed to help with this, but we'll see.


> I've written code to do so many times -- Numerical Recipes provides the
> basic idea in its matrix-packing routines.  Simply allocate e.g.
> 
> double **..*m,*v;

I understand the mechanics of it. I'm not sure what my concern was;
I thought about this a few months ago and decided there was some
semantic problem. But it looks like the **..*m has exactly the
same life-cycle and maintenance as *v, so I'm not sure what I
was thinking. It might come back to me.


> I can probably dig this code out and post it here (or post a link to it)
> if you want to look at it.  I think I went up to 9 or 10 dimensions (way
> more than I needed) and did a lot of things via brute force to make the
> purely mechanical recursion obvious; I planned an eventual rewrite into
> a more compact form that I never got around to.

Sure. Feel free to post it. I like to look at everything and
cherry-pick the best parts from other people's stuff. Makes
things a lot easier.


--
G. Jungman

Re: gsl container designs

[Brian Gough]

At Wed, 06 Jan 2010 13:44:56 +0200,
Tuomo Keskitalo wrote:
> Do you mean the compatibility to GSL 1 types by the compatibility cost? 
> When talking about GSL 2 I don't think we should give too much value to 
> maintaining backwards compatibility. GSL 1 is not going to cease to 
> exist, and people who have tied themselves deeply to GSL 1 data 
> structures can continue to use it.

I'd say it's a significant benefit if people don't have to modify
their code. 

-- 
Brian Gough

Re: gsl container designs

[Tuomo Keskitalo]

Hello,

On 11/25/2009 02:54 AM, Gerard Jungman wrote:

> Here are header files for a couple different approaches to containers.
> I didn't bother with any implementations; it seems obvious how to
> implement most of these functions.
> 
> The designs are not complete, but they express most of
> the important stuff.

Your design 1 looks very interesting. I don't think your designs would 
support banded and packed matrices, or would it? (see the last pages of 
"Legacy BLAS: C Interface to the Legacy BLAS" at 
http://www.netlib.org/blas/blast-forum/)

I could not think of anything why your design 1 would not work for me. 
The index mapping is particularly elegant!

-- 
Tuomo.Keskitalo@iki.fi
http://iki.fi/tuomo.keskitalo

Re: gsl container designs

[Gerard Jungman]

On Wed, 2010-01-06 at 14:04 +0200, Tuomo Keskitalo wrote:

> Your design 1 looks very interesting. I don't think your designs would 
> support banded and packed matrices, or would it? (see the last pages of 
> "Legacy BLAS: C Interface to the Legacy BLAS" at 
> http://www.netlib.org/blas/blast-forum/)

Thanks. It's been a long time since I looked at that document.

You are right; I have only considered the dense case,
purely from laziness. However, I think the banded and
packed cases can be handled similarly. One needs only
to define the indexing schemes to layer on top of
the (essentially universal) data structs. Finally,
add a syntactic layer, to inject appropriate names,
analogous to the syntactic layer that declares
vectors, matrices, etc.

At least, I think it is not conceptually difficult.
Maybe I am missing something.

--
G. Jungman