This is the mail archive of the gsl-discuss@sources.redhat.com mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]

Re: gsl Development Query


Gerard Jungman wrote:

> Randall Judd wrote:
> 
> > Part of the specification is the interface
> > which allows people
> > to bring data into or export data from VSIPL.
> 
> This is an important point.
> It may be that data interchange is issue number
> one for this kind of library.  Certainly, I feel that
> it is the number one issue for vector/matrix/etc.
> implementations in numerical libraries.
> No library can afford to act as an island fortress.
> As a general comment, I think that generic programming
> (which for all practical purposes
>  means algorithm templates in C++)
> offers the best hope at this point for freeing ourselves
> from the details of data representations, in the wide sense.
> I see no real solution,
> when restricted to the context of a C library like GSL;
> we just have to find the right compromise,
> given the right set of restrictions on the problem domain.
> 
> By the way, I looked at the VSIPL draft
> to see what you were doing about this,
> but I wasn't sure what the appropriate section was.
> From what I did see,
> I guessed that your block/view separation was the key
> to data interchange, allowing conformant views
> to be layered over imported data blocks.
> Is this what you meant,
> or did you have something more specific in mind?
> 
> I certainly agree with this approach.
> GSL actually works this way as well
> (contrary to whatever Tisdale seems to be going on about),
> though it is not as complete.
> Certainly I hope that by the time we are done
> we will have a similar level of completeness,
> addressing such issues as views
> layered over different underlying complex array types.
> It's not hard to do, and it fits in with what we have so far.
> Brian, are you listening?

The VSIPL API was designed to accommodate
embedded Digital Signal Processors (DSPs).
These machines typically have memory spaces
which can't be accessed through a standard ANSI C pointer.
The memory may be in vector registers, cache,
program memory space, data memory space,
other processing nodes in a multiprocessor system, etc.
A VSIPL block object allows VSIPL functions
to reference these non standard memory spaces
as well as the memory space that can be referenced
through an ordinary ANSI C pointer.
Apparently, GSL blocks can only reference memory
that can be referenced through an ordinary ANSI C pointer.
You can always obtain an ordinary ANSI C pointer
to the underlying data array from a GSL block or view object
but you can't always get such a pointer
from a VSIPL block object.

I assume that you expect all of the GSL users
are like the users at LANL who mainly program
super computers and high performance workstations.
If so, you will probably get away with the "flat"
virtual memory model that the GSL currently supports.
But if you think that GSL users might want to port
applications to some of the new high performance DSPs,
you might want to reconsider the VSIPL architecture.

> Of course, there are many hidden problems.
> Especially that, in principle,
> one needs to optimize algorithms differently
> if using a certain type of view
> implies a poor strategy for memory accesses.
> This rapidly becomes complicated, parametrizing
> both the data layouts and algorithms together.
> But that makes it interesting too.
> 
> > The main problem with FFTW is that
> > they use an array of complex.
> > Most of the time an array of complex
> > looks like interleaved complex
> > (which is, I think, what GSL uses for complex arrays).
> > For the user interface,
> > VSIPL supports both interleaved or split arrays
> > (internally a complex block is just a block of complex).
> 
> Just a minor technical note,
> for readers who may not be aware of the details.
> As I'm sure you know, this layout of complex numbers
> (coincident values of underlying type, like double x[2])
> is standard in fortran and (I believe) the C99 draft.
> From a practical standpoint, I think there is no other way to go,
> because you have to be able to inter-operate with old code,
> especially old fortran code.
> The fact that C++ does not specify the implementation details
> for std::complex is, I believe, a deficiency of the C++ standard.
> Eventually they will have to come around and fix this;
> we can't glide along forever, relying on the the fact that
> most compiler/platform combinations end up with
> 
> 	struct { double x; double y; }
> 
> and
> 
> 	struct { double x[2]; }
> 
> aligned "by accident" on the same boundaries.
> 
> If anybody wants to clarify
> what is going on with this issue out there,
> I would like to know myself.
> This was one of the stickier points
> that Brian and I had to make a decision on.

Application programmers who aren't particularly concerned
about the portability of their code needn't worry about this.
Actually, virtually every ANSI C compiler
will store both of these types in exactly the same way.
The problem is that the ANSI C standard does not guarantee
that they will be stored in exactly the same way
so the second type is safer.

The C++ standard does not distinguish
Cartesian complex numbers from polar complex numbers
or any other representation of complex numbers.
High performance numerical computing, on the other hand, has always
used Cartesian complex numbers represented by a pair of real numbers
in an array of two elements with the imaginary part following
immediately after the real part (with unit stride between them.)
This representation of complex numbers predates Fortran
as does the notion of real and complex multidimensional views
of a one dimensional array of real numbers.
These notions of views were incorporated into the design
of the Fortran programming language using the EQUIVALENCE declaration.

> > I guess my point is that you probably should have your own FFT's.
> 
> Yes, I can't disagree.
> If somebody wants to drop multi-dim FFTs in our lap,
> that would be great. I just don't want to have to do it.
> And if we take source code from somewhere else, we have a choice.
> We can figure out how to create an appropriate layer
> that insulates us, or we can apply the old copy/paste methodology.
> Now, I find copy/paste to be the most revolting form of reuse,
> as I'm sure most people do.
> On the other hand, it all boils down to two questions:
> 
>   (1) Is the source you are snarfing stable enough?
>   (2) Can you imagine ever wanting to swap something else in?
> 
> In the cases at hand, the answer to (1) is probably yes.
> But I have never yet found a case where the answer to (2) was no.
> That's somewhat subjective, but that's how I feel about it.
> 
> Now, Brian may object to the copy/paste statement. After all,
> what one of us would really do would be to read about it,
> look at a couple different implementations,
> and then code something up as cleanly as possible.
> But in my book this is just a more sophisticated
> (and not necessarily more intelligent) form of copy/paste.
> 
> > It may not be reasonable to make an interface to FFTW
> > and keep the other interfaces in your library coherent.
> 
> My feeling is that, in this case, it won't be hard.
> Similarly, it would not be hard to use
> an underlying fortran implementation for the BLAS library.
> All you have to do
> is specify an appropriate middle layer in each case.

Maybe.  Maybe not.  You can never be sure that
you will have compatible Fortran and C compilers
on any given platform.
You might need to write an interface in assembler
to call Fortran subroutines from a C program.

> What I can say for sure is that
> if we create an appropriate middle layer first,
> then we have options, and if we don't
> then it could end up being very difficult
> for us or anyone else to put something new in later.
> 
> In this regard, this is a question about modularity
> with respect to the implementations.
> If you can make things more flexible
> without sacrificing simplicity
> at the highest level of abstraction
> (where most users will live),
> and it doesn't cost too much to do so, then why not?
> 
> > If somebody really needs FFTW speed
> > and also wants to use your library,
> > then they should figure out how to make them interface.
> > Hopefully, you give them enough information so they can do this.
> 
> Well, maybe we can kill all these birds
> by making the appropriate middle layer work ourselves.
> It can help us, because we get to link in
> a working implementation with minimal effort,
> and it can help other people if it saves them the effort.

I don't think that you need to get bogged down
in the GSL implementation details.
It may seem like a huge project to you now
but eventually you will get some perspective
and realize that they really are just details.
A portable version of the GSL is certainly a good idea
but there will never be
just one really good implementation of the GSL.
If the GSL is to survive, then, just like the VSIPL,
there must be several implementations each tailored
to a specific platform or class of computer architectures
and/or application domains.

Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]