Re: gsl Development Query

Re: gsl Development Query

[E. Robert Tisdale]

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.

Re: gsl Development Query

[Eleftherios Gkioulekas]

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

You may check that with Autoconf. At worst, you can fall into f2c as
a last-resort backup. 

-lf.

Re: gsl Development Query

[Gerard Jungman]

"E. Robert Tisdale" wrote:
> 
> Apparently, GSL blocks can only reference memory
> that can be referenced through an ordinary ANSI C pointer.

Yes. That was an assumption. GSL was never targetted
to more esoteric platforms.


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

Or how about just non-descript intel boxes,
like the one I am sitting in front of? The fact
is that the audience for GSL is precisely that
(very large) group of people who sit in front
of boxes just like the one in front of me.


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

Well, "particularly concerned" has different meanings
for different people. I think being insulated from
issues like the alignment is important. People learned
that as soon as 64 bit boxes became common.


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

Uhh... exactly how many lines of code have you
contributed to GSL lately, that you are in
a position to measure anything related to
my effort?


-- 
G. Jungman

Re: gsl Development Query

[E. Robert Tisdale]

Steve ROBBINS wrote:

> OK, then why do you care to improve the API?  
> 
> Are you using GSL or do you plan to?
> That seems unlikely,
> as your previous message stated that
> you have already "done GSL right". 
> 
> I'm just wondering what your interest in GSL is.

I think that I've said this before.
I'm interested in a good general purpose numerical library
just like all of the other contributors to this discussion.
I had hoped that the GSL would be that library.  It isn't.
And it may be too late to get it turned around
and headed in the right direction.  I don't know.

Re: gsl Development Query

[E. Robert Tisdale]

Brian Gough wrote:

> Randall Judd writes:
> 
> > In VSIPL every view contains a pointer to the block. 
> 
> I think we can accomodate this in GSL
> by not setting the block field to NULL for views
> and having a separate field in the struct to indicate ownership
> of the memory instead.  Then one would always be able
> to determine the underlying block for any view.
> 
> This seems straightforward to support, so I will do that --
> it only needs a change in the struct definition
> and a few associated routines.

Do you really need to tell us these details?
It's your implementation.  Do what you want.
But give us library functions that return

  1.) a pointer to the block,
  2.) the offset from the beginning of the block
      to the first element of the vector,
  3.) the length (extent) of the vector and
  4.) the stride between elements of the vector.

You can implement them as inline functions
or even C preprocessor macros if you like.
Take a look at
The ANSI C Numerical Class Library

	http://www.netwood.net/~edwin/svmt/

and read cncl/src/vector/vector.hP for an example
of how you might implement them.

Re: gsl Development Query

[E. Robert Tisdale]

Steve ROBBINS wrote:

> ... so remind us again why you subscribed to gsl-discuss
> and proceeded to berate the developers of GSL?
> 
> -S
> 
>         "you catch more flies with honey than with vinegar"

I do not berate the GSL developers.
I berate the GSL API.
Professional people must learn
not to take criticism of their work personally.
My work (and Randy's work) on the VSIPL API standard
received a great deal of criticism
but I think it helped me and it made the VSIPL API better.
I only hope that the GSL API will benefit as much
from the criticism that I offer.
Just try to think of it as peer review.

Re: gsl Development Query

[Steve ROBBINS]

On Mon, 4 Sep 2000, E. Robert Tisdale wrote:

> Steve ROBBINS wrote:
> 
> > ... so remind us again why you subscribed to gsl-discuss
> > and proceeded to berate the developers of GSL?
> > 
> > -S
> > 
> >         "you catch more flies with honey than with vinegar"
> 
> I do not berate the GSL developers.
> I berate the GSL API.

You sidestep the question.  Why do you berate the API?

> I only hope that the GSL API will benefit as much
> from the criticism that I offer.

OK, then why do you care to improve the API?  

Are you using GSL or do you plan to?  That seems unlikely, as your
previous message stated that you have already "done GSL right". 

I'm just wondering what your interest in GSL is.
-S

Re: gsl Development Query

[E. Robert Tisdale]

Randall Judd wrote:

> I don't attach VSIPL blocks to GSL views.
> The idea is to associate (a word I prefer to attach)
> VSIPL blocks with GSL blocks (not views).
> Then define VSIPL vector and/or matrix views
> which map the block data space
> the same way that the GSL views map the GSL blocks data space.

A VSIPL block is like a C++ valarray.
The are both pure abstractions.
You can read what Kent Budge (kgbudge@valinor.sandia.gov)
has to say about valarray in his article
"Re: OONSTD: to early for standards?" at

	http://www.oonumerics.org/oon/oonstd/archive/0018.html

Both the VSIPL block and the C++ valarray
give the application programmer the "illusion"
that the underlying data array
is a contiguous block of real numbers
regardless of how they are actually stored in memory.
The C++ compiler may actually conspire
with the operating system and the computer architecture
to create this illusion for application programmers
who use valarray and return a "pointer" to an element
in a valarray
The ANSI C compiler is not expected to participate
in this illusion for VSIPL block objects.
The actual data representation of VSIPL block objects
is hidden from the application programmer
and, in general, elements can only be accessed
through VSIPL functions.

A GSL block isn't abstract at all.
It is defined to be

        struct gsl_block_struct {
          size_t size;
          double* data;
          };
        
        typedef struct gsl_block_struct gsl_block;

Aplication programmers are permitted to access
the double* data member directly and they are guaranteed that
they can always obtain a pointer to the underlying data array.
This means that the GSL can never distribute the data array
among the nodes of a multiprocessor system for example.

> As I understand it, your GSL block contains some data space.
> Then you slice this data space up.
> The first of your vector views creates
> the block, data space, and a vector.

All this does is to combine two VSIPL operations.
The first vector object "owns" the block.
The block should never be destroyed
until just before the owner is destroyed.

> The vector view contains a pointer to the block
> plus a pointer to the beginning of the data.
> When you take of subview of the vector,
> you set the block pointer to NULL,

This is just a precaution to help ensure
that the block is never destroyed
until just before the owner is destroyed.

> get a pointer to the beginning of the subviews data space

Essentially, they just precompute vector.data = block.data + offset
They could obtain the offset = vector.data - block.data
if they hadn't discarded the pointer to the block.
It's too bad really.  If they had retained the pointer to the block,
they could have used it to help determine
whether source and destination views overlap or not.

> and then set the attributes of the subview
> to the proper stride and length for the new vector
> which must not exceed the data space of the block.
> 
> In VSIPL every view contains a pointer to the block.
> The block contains the entire possible data space.
> For a GSL vector view attached to a block,
> the view contains a pointer to the beginning of data,
> a length equal to the data space length,
> and a stride equal to one.
> The corresponding VSIPL vector view contains
> an offset of zero, a stride of one
> and a length equal to the length of the block.
> 
> So, when I create a block,
> I would want the block to encompass the entire GSL block.
> 
> Now I could create a new block to do the GSL subview,
> but the new block would encompass some of the same data space
> as the previous VSIPL block.
> This would be an error from VSIPL's point of view.
> What I do in VSIPL is set an offset into the data space
> so my subview points to the same data space as the GSL subview.
> Then the stride and length would be the same as GSL.

Re: gsl Development Query

[E. Robert Tisdale]

Randall Judd wrote:

> The Core, and Core Lite documents are what are called profiles.
> No vendor, except perhaps SKY, has completed all of VSIPL yet.
> But they wanted to be able to say they are VSIPL compliant.
> So we defined profiles to give them a starting point
> and a way to claim compliance to something.
> Core Lite is only about 129 function
> and is a very easy starting point for most vendors.
> A great many, surprisingly so, current Navy applications
> appear to be mappable to just a core lite implementation.
> The core profile has about 500 or so functions
> and include more complicated things.
> In the future, there may be more profiles.
> A profile is basically a defined set of functions
> that has been approved by the VSIPL forum.

It is important to distinguish between
what a standard specifies and what it mandates.
The problem with the BLAS library API standard, for example,
is that it mandates everything that it specifies.
You can't implement part of the BLAS library API standard and
say that you have implemented the BLAS library API standard.
Somebody would get upset and accuse you of fraud.
The BLAS library API standard designers couldn't get
BLAS library developers to implement everything
that they wanted in the BLAS library API standard
so they compromised.  They left out a lot of functionality.
The result was a much impoverished API design
from the numerical application programmers point of view.
You can check this out by reading
"A Set of Level 3 Basic Linear Algebra Subprograms,"
by Jack Dongarra, Jeremy Croz, Iain Duff and Sven Hammarling,
Section 8 Rationale which you can download
as a PostScript file blas3-paper.ps from

	http://www.netlib.org/blas/

The VSIPL profiles decouple specifications from mandates.
They "partition" the VSIPL library into manageable chunks
which depend upon each other in a natural way.
Matrix objects depend upon vector objects,
linear algebra depends upon matrix objects,
tensor objects depend upon matrix objects, etc.
so library developers can build a complete VSIPL library
in increments over time.

Still, the VSIPL API standard is limited in scope.
The vector part includes some linear algebra
but it is mainly intended for digital signal processing.
Eventually, it is supposed to include
digital image processing as well.
It was never really intended to be
a general purpose Scientific Computing library
as I (and Randy I suppose) hoped the GSL was supposed to be.

Re: gsl Development Query

[Randall Judd]

I guess I don't know if Bobs comments belongs on this discussion group
since it does not talk very much about GSL but their are a couple of
corrections I would like to make.

First a comment. GSL keeps everything out front. VSIPL hides everything. If
I was doing GSL I very well might write my own support functions. By that I
mean I don't have to make a GSL block using a GSL funciton. I could just do
gsl_block a_block;
double data[100]
a_block.data = data;
a_block.size = 100;

You can not do this in VSIPL. So part of the job of a profile is to make
sure all the correct support functions are included to handle all the
functions you want to do. 

>
>The VSIPL profiles decouple specifications from mandates.
>They "partition" the VSIPL library into manageable chunks
>which depend upon each other in a natural way.
>Matrix objects depend upon vector objects,
This is wrong. Matrix objects stand alone and do not depend on vector objects.

>linear algebra depends upon matrix objects,

Dito for tensors.
>tensor objects depend upon matrix objects, etc.
>so library developers can build a complete VSIPL library
>in increments over time.
>

Don't know what the following means. We defined some functions. We left a
lot out. We still have a lot of functions defined. 

I think one of the main differences between VSIPL and GSL is it is easier
to add new functionality to GSL. If somebody gives you a well written
function following the rules as set out by your forum you would probably
just add it in. I expect you are not so much concerned with adding new
functions, it's writing them to start with. With VSIPL it is convincing
vendors that they should all write this function so they are not likely to
put one in just because somebody has written it already. To get new
functionality into VSIPL requires a lot of work, and I don't mean programing.

>Still, the VSIPL API standard is limited in scope.
>The vector part includes some linear algebra
>but it is mainly intended for digital signal processing.
>Eventually, it is supposed to include
>digital image processing as well.
>It was never really intended to be
>a general purpose Scientific Computing library
>as I (and Randy I suppose) hoped the GSL was supposed to be.

I actually don't have any hopes for GSL except I hope you succeed in your
efforts. It sounds like a lot of work, and I know because I still have a
lot of work to do on my VSIPL implementation. My main concern for GSL is
some API inconsistency. For instance shouldn't

gsl_fft_complex_radix2_forward(gsl_complex_packed_array *data,size_t
stride, size_t n)

have some sort of vector input? I am not sure why you did it this way. Note
that I don't really know all that much about your library. But I get the
feeling that every programer is defining his own interface. You may end up
with a lot of nice functions but they won't fit together very well and will
be hard to program with.

                      Randy Judd

Re: gsl Development Query

[Brian Gough]

 > My main concern for GSL is some API inconsistency. For instance shouldn't
 > 
 > gsl_fft_complex_radix2_forward(gsl_complex_packed_array *data,size_t
 > stride, size_t n)
 > 
 > have some sort of vector input?

This is one of several things I need to document, I will try to do
that. There are some variations in the API are part of the design and
there are some inconsistencies which are genuine, due to different
developers following different paths and which we are still working
on.  I will explain this particular point here since it is actually
part of the design,

For purely "1d" functions we use the C-style arguments (double *,
stride, size) so that it is simpler to use the functions in a normal C
program, without needing to invoke all the gsl_vector machinery.

The philosophy here is to minimise the learning curve. If someone only
needs to use one function, like an fft, they can do so without having
to learn about gsl_vector.  

This leads to the question of why we don't do the same for matrices.
In that case the argument list gets too long and confusing, with
(size1, size2, tda) for each matrix and potential ambiguities over row
vs column ordering. In this case, it makes sense to use gsl_vector and
gsl_matrix, which take care of this for the user.

So really the library has two levels -- a lower level based on C types
for 1d operations, and a higher level based on gsl_matrix and
gsl_vector for general linear algebra.

Of course, it would be possible to define a vector version of the
lower level functions too. So far we have not done that because it was
not essential -- it could be done but it is easy enough to get by
using the C arguments, by typing v->data, v->stride, v->size instead.
A gsl_vector version of low-level functions would mainly be a
convenience.

Re: gsl Development Query

[E. Robert Tisdale]

Gerard Jungman wrote:

> Uhh... exactly how many lines of code
> have you contributed to GSL lately,
> that you are in a position to measure anything
> related to my effort?

Well, now that you brought it up,
exactly where did you get the idea
for the GSL block, vector view and matrix view objects.
I first proposed the VSIPL block, vector view,
matrix view and tensor view objects
in a presentation before the VSIP Forum in April, 1997.
I implemented the original VSIP portable reference library
which served as the starting point for the implementation
by Randy Judd which you can get on the internet today

    http://www.vsipl.org/

I am the author of
The C++ Matrix class library

    http://www.netwood.net/~edwin/Matrix/

which has been available free of charge on the internet
for over a decade now.  Hundreds of people use it.

I am the author of
The C++ Scalar, Vector, Matrix and Tensor class library

    http://www.netwood.net/~edwin/svmt/

which has been available free of charge in the internet
for over two years now.   Hundreds of people use it.

I am the author of
The ANSI C Numerical Class Library

    http://www.netwood.net/~edwin/svmt/

which is basically the GSL done correctly.

I've contributed.
I've written many more lines of code than you have
and distributed them free of charge for well over a decade now.
If you choose to ignore them, that's your problem not mine.

It isn't any wonder that you have trouble
getting people to contribute to the GSL.
You have no plan.  You have no design.  You have no goals.
At least none that you have or can articulate.
I have very little confidence in the GSL right now.
I don't think that it is going anywhere.
And I don't think anyone else does either.
I don't want to waste my time on a project
that is doomed from the start
and I don't think that anyone else does either.

Please, if you are serious about the GSL,
go back and read my earlier posts to this discussion
and try to answer some of the questions that I asked.
It really will help.

Re: gsl Development Query

[Steve ROBBINS]

On Sat, 2 Sep 2000, E. Robert Tisdale wrote:

> I am the author of
> The C++ Matrix class library [...]
 
> I am the author of
> The C++ Scalar, Vector, Matrix and Tensor class library [...]
 
> I am the author of
> The ANSI C Numerical Class Library
> 
>     http://www.netwood.net/~edwin/svmt/
> 
> which is basically the GSL done correctly.

... so remind us again why you subscribed to gsl-discuss and proceeded to
berate the developers of GSL?

-S

	"you catch more flies with honey than with vinegar"

Re: gsl Development Query

[E. Robert Tisdale]

Eleftherios Gkioulekas wrote:

> At worst, you can fall back onto f2c as a last-resort backup.

Provided that you have the Fortran source code.
You might not be able to use the proprietary high performance version
of the BLAS library that is available on your platform.
You might be obliged to fall back onto the BLAS portable reference library
and take an unacceptable performance hit.

Re: gsl Development Query

[E. Robert Tisdale]

Gerard Jungman wrote:

"E. Robert Tisdale" wrote:

> Actually, counting by man years of effort,
> most of the "programmers at LANL"
> spend their lives writing very expensive,
> large, and nondisposable applications.
> 
> It's fairly rude to tell other people what they do.

I'm familiar with LANL.
Most of the programmers at LANL are scientists and engineers.
They don't think of themselves as programmers.
They don't get paid to write programs
but they do write programs to help them do their own work.
They seldom write programs that other people use.
Usually, their programs are of no use once they have completed
the task that their programs were written to perform.
Their programs are usually archived for a while
and then deleted when nobody can remember what they do
and the storage space is needed for something else.

> Apparently you consider yourself a "professional programmer",
whatever that is.

Someone who gets paid
specifically to write programs for other people to use.

> Ok, why don't you put a little substance behind your mouth
> and contribute some "professional code" to this project?
> Maybe you could earn the right to badmouth the people
> who have actually done something.

Take a look at
The C++ Scalar, Vector, Matrix and Tensor class library

	http://www.netwood.net/~edwin/svmt/

Take a look at
The ANSI C Numerical Class Library on the same page.

Take a look at
The Vector, Signal and Image Processing Library

	http://www.vsipl.org/

I've contributed.

> I think most "professional programmers"
> would fail to find their own assholes in the dark.
> I can't be expected to account for their behavior.
> 
> The question at issue here is not one of API design.
> The main point of this discussion
> is how to most efficiently implement the functionality
> that is required to reach the stated goals of the project.
> Don't change the subject, and don't insult the people
> who are actually doing something.

Personally,
I don't see anything wrong with Brian's implementation.
It's his implementation and he should be left alone
to develop it as he sees fit.
If you don't like it, you should re-implement it yourself.
The problem is that Brian isn't simply implementing
a standard API like the BLAS library API.
He is inventing a new API and
his API is designed around his implementation.
Eventually, he will discover that he has made mistakes.
Everybody makes mistakes but if application programmers
start using the GSL, it will soon be practically impossible
to correct those mistakes without breaking
all of the existing application programs
that depend upon the GSL.

> Finally, I apologize to the readers of this list
> for even bothering to reply to these messages.
> Perhaps doing so is simply a waste of everybody's time.
> If somebody out there has something intelligent to say
> about the implementation policies for this project,
> please do not be discouraged by this exchange from doing so.

If you invite discussion in an open Forum like this,
you are bound to read things that you don't like.
Right now the GSL API is pretty awful.
The so-called error handling is ineffective
and the API provides many more opportunities
to make programming errors than any error handling scheme
could ever hope to deal with.
It is difficult and dangerous to use.
It adds unnecessary overhead to numerical functions.
There doesn't appear to be any design rationale
behind the GSL that anyone can articulate.
In short, the GSL API stinks.
Numerical application programmers should avoid the GSL
until the GSL designers get their act together
and publish a document which specifies the GSL API
and the rationale behind their design decisions.

Re: gsl Development Query

[E. Robert Tisdale]

Gerard Jungman wrote:

> Look, we can all sit and play C++ compiler.  But there's no point.
> The implementation language for this project does not make a distinction
> and nobody at this end is going to make pedantic distinctions.
> People who use C libraries have to apply some common sense.
> That is a fact (unfortunate perhaps) of the language.
>
> And he did it because it was useful to write it down somewhere.
> Don't quibble about where he wrote it down.
>
> Nobody is painting any corners
> and nobody is obliged to do anything.

You are wrong.
The documentation specifies each of the data members in

        typedef struct {
          size_t n;
          size_t nf;
          size_t factor[64];
          gsl_complex *twiddle[64];
          gsl_complex *trig;
          double *scratch;
          } gsl_fft_wavetable_complex;

Application programmers can reference and modify
the private data members
of a gsl_fft_wavetable_complex object directly.
Apparently, the documentation guarantees
that a gsl_fft_wavetable_complex object
will always have this data representation and
the data members will always have this interpretation.
That means that the GSL developers
can never change the data representation
or the implementation of the GSL FFT
without breaking application programs
that depend upon this representation
and the implementation that it implies.

The GSL provides no library functions
which application programmers can use
to access the private data member indirectly.
The GSL makes no attempt to hide
the actual data representation
from the application program.
The GSL documentation makes no distinction
between the API and implementation details.
Application programmers can't "apply common sense"
because GSL application programmers
and GSL library developers don't share
the necessary information.
Application programmers would need to be clairvoyant!
And upon reading your mind and Brian's mind,
they would discover that you and he can't even agree
so how would you expect application programmers to decide
what is part of the API and what is implementation detail?

I'm not talking about anything exotic here.
All that I am talking about is simply applying
good programming practice.
It has nothing to do with any particular programming language.
Pick up a copy of "Data Structures and Algorithms"
by Alfred V. Aho, Jeffrey D. Ullman and John E. Hopcroft
or just ask any sophomore computer science student
who has recently completed a course in data structures.

Re: gsl Development Query

[E. Robert Tisdale]

Gerard Jungman wrote:

> Brian Gough wrote:
> 
> > What you've outlined is something much bigger
> > than GSL (i.e. it's not GSL).
> 
> No. What I have described is an efficient way
> to obtain the required functionality in GSL.
> These are implementation details.
> If you feel there are some fundamental flaws
> in existing implementations,
> then you are free to reimplement. I don't see
> the flaws in existing (and Free!) implementations
> such as BLAS and FFTW.
> 
> And how can it possibly matter to the end user?
> All they want is something that works,
> and that works in the right way,
> without making undesirable assumptions
> about usage patterns or data representation.
> These important issues are _interface_ issues.
> Their only relation to the elemental entity
> turning the crank under the hood
> is to specify certain constraints
> that any such implementation must satisfy.
> If you find a solution to those constraints,
> then you have satisfied the design criteria.
> If the solution happens to already exist,
> then so much the better; you can go home early.
> 
> > GSL is just a collection of routines for numerical computing,
> > written in C with modern coding conventions.
> 
> Yes, that's what it says in The Big Book of GSL.
> Now, rather than quoting the precious scripture to me,
> why don't you apply your native intelligence
> and think about whether or not it makes any sense.
> 
> The so-called "design" of this library
> was never anything more than some vague wish list.
> The original idea that you could provide
> "loosely coupled" libraries for the different tasks
> was wrong.  This was made clear
> when Brian and I got down to worrying about
> the various multi-dimensional functionalities
> which depend explicitly on data structures
> for vectors and matrices and linear algebra support.
> The linear algebra bottleneck should be obvious to everyone
> who has done any numerical computing;
> if you don't solve this problem first, you are screwed.
> 
> When you take the time to look at the requirements
> for implementing higher level algorithms,
> you realize that the design must be layered
> and you must solve the fundamental problems first.
> These include things like control of and access
> to the floating point capabilities of the machine,
> data representation issues, and basic linear algebra.
> 
> When you (Brian) and I faced these problems squarely,
> about two years ago, we found it necessary
> to redesign and reimplement large sections of the project,
> fixing things related to the issues above,
> as well as doing the right thing at higher levels,
> such as introducing a consistent framework for iteration.
> 
> When you say "well, we're just a collection of tools",
> you are just echoing the original design goals,
> which have been demonstrated, by your own work,
> to be confused.
> 
> > Maybe a "universal interface/library"
> > of the type you suggest is worthwhile --
> > but it would be a different project from GSL.
> 
> I'm not talking about a "universal interface/library",
> whatever that is.  I'm talking about the right way
> to provide an underlying implementation
> for some of the required functionality.
> 
> > If you need wrappers for fortran codes,
> > they are available in other places
> > (e.g. Octave's library, NumPy, ...).
> 
> Yes, they are. So?
> 
> Nobody "needs" wrappers for fortran codes.
> They need tools that do certain tasks.
> Wrapping fortran just happens to be
> a convenient way to get certain functionality.
> That's why people do it,
> not because they like to know that
> a little fortran monkey is turning the crank
> in their greasy little simulated worlds.
> Who cares what turns the crank, as long as it turns?
> 
> > The idea of GSL is that wrapping
> > is never completely satisfactory
> > and so it's better to rewrite things.
> 
> Let me paraphrase this for the general reader:
> The idea of GSL is to write half-assed implementations
> of some ill-defined collection of numerical foo,
> sometimes of functionality that already exists,
> with no real thought going into the overall design,
> and to pat each other on the back
> for having done everything the hard way
> and put it under the GPL.
> Well, GPL'ed foo is still foo;
> that's one thing I know for sure
> after having watched the free software world work
> for the last ten years.
> 
> You have this idea that GSL is just some thing
> that we are tinkering with in the backyard,
> not really "important" enough to do the right way;
> after all, nobody's going to the moon riding on a GSL rocket.
> This is insulting to everyone involved, including yourself.
> 
> > However much you wrap unreadable fortran,
> > it's still unreadable fortran.
> 
> Huh? What does this mean? This is a non sequitur.

Evidently, you and Brian
each have a very different vision for the GSL.
Brian's aspirations are rather modest.
He just wants to implement a numerical library
that he can distribute free of charge to people
who can't afford the $50 or whatever it costs
to purchase a copy of Numerical Recipes in C.
For him, the GSL is just another software project.
His plan is to simply allow the GSL to evolve.
Apparently, he believes that GSL users
are amateur application programmers who want to write
short lived disposable application programs
like most of the programmers at LANL for example.
He doesn't worry about supporting poor API design decisions
because he doesn't expect any applications using the GSL
to survive much beyond the next version of the GSL.

Personally, I would prefer an API for the GSL that
professional as well as amateur programmers could use.
Professional programmers are more interested in writing
long lived reusable application programs.
I don't think that many professional programmers
would dare risk using the GSL because, apparently,
it isn't reasonable to expect Brian or anyone else
to support it for very long.

Anyway, I don't think Brian is prepared to design an API.
If you want someone to design an API for the GSL,
you will probably need to do it yourself which means that
you will need to contribute to the design documentation
and submit specific proposals -- in writing.
The API should not dictate how it is to be implemented.
Implementation details should be left up to
the library developer --
Brian if he agrees to implement the API.

Re: gsl Development Query

[Gerard Jungman]

"E. Robert Tisdale" wrote:
> 
> Evidently, you and Brian
> each have a very different vision for the GSL.
> [...]
> Apparently, he believes that GSL users
> are amateur application programmers who want to write
> short lived disposable application programs
> like most of the programmers at LANL for example.

Actually, counting by man years of effort, most of
the "programmers at LANL" spend their lives writing
very expensive, large, and nondisposable applications.

It's fairly rude to tell other people what they do.


> Personally, I would prefer an API for the GSL that
> professional as well as amateur programmers could use.
> Professional programmers are more interested in writing
> long lived reusable application programs.

Apparently you consider yourself a "professional programmer",
whatever that is. Ok, why don't you put a little substance
behind your mouth and contribute some "professional code"
to this project? Maybe you could earn the right to
badmouth the people who have actually done something.


> I don't think that many professional programmers
> would dare risk using the GSL because, apparently,
> it isn't reasonable to expect Brian or anyone else
> to support it for very long.

I think most "professional programmers" would fail
to find their own assholes in the dark. I can't
be expected to account for their behaviour.


> Anyway, I don't think Brian is prepared to design an API.

The question at issue here is not one of API design.
The main point of this discussion is how to most efficiently
implement the functionality that is required to reach
the stated goals of the project. Don't change the
subject, and don't insult the people
who are actually doing something.


Finally, I apologize to the readers of this list
for even bothering to reply to these messages.
Perhaps doing so is simply a waste of everybody's
time. If somebody out there has something
intelligent to say about the implementation
policies for this project, please do not be
discouraged by this exchange from doing so.


-- 
G. Jungman

Re: gsl Development Query

[E. Robert Tisdale]

Brian Gough writes:

> E. Robert Tisdale writes:
> 
> > The GSL should allow application programmers
> > to create FFT objects like the "plan" object in FFTW.
> > The overhead of creating an FFT object is amortized
> > over several invocations of an FFT of a particular length.
> 
> Yes, this is how it works --
> see the documentation for gsl_fft_complex_alloc
> in the manual for details.

The FFT object (what you call a gsl_fft_wavetable_complex) is OK.
But then you go on to itemize the private members:

	typedef struct {
	  size_t n;
	  size_t nf;
	  size_t factor[64];
	  gsl_complex *twiddle[64];
	  gsl_complex *trig;
	  double *scratch;
	  } gsl_fft_wavetable_complex;


Why did you do this?
It isn't necessary for the application programmer
to have any of this information.
You present it in the manual as if it is part of the API.
Did you really need to make a commitment
to this particular algorithm
and this particular data representation?
I don't think so.
I think that you've painted yourself into a corner
and you will be obliged to support this particular
data representation and the implied algorithm forever.

There is almost never any one best algorithm
much less one best implementation of anything
in numerical computing.
No matter how gifted and careful you are,
your implementation and/or data representation
is likely to be suboptimal on some platform.
It is best to hide as many details
of the data representation and/or implementation
so that you or some other programmer can re-implement
the functionality without breaking
any of the applications which depend upon it.

If, on the other hand, you only meant to describe
the data representation and algorithm
that you were implementing, you should have
separated it from the specification of the API -- 
preferably, in another document.

Re: gsl Development Query

[Gerard Jungman]

"E. Robert Tisdale" wrote:
> 
> The FFT object (what you call a gsl_fft_wavetable_complex) is OK.
> But then you go on to itemize the private members:
> 
>         typedef struct {
>           size_t n;
>           size_t nf;
>           size_t factor[64];
>           gsl_complex *twiddle[64];
>           gsl_complex *trig;
>           double *scratch;
>           } gsl_fft_wavetable_complex;
> 
> Why did you do this?
> It isn't necessary for the application programmer
> to have any of this information.

Look, we can all sit and play C++ compiler. But there's
no point. The implementation language for this project
does not make a distinction, and nobody at this end is
going to make pedantic distinctions. People who use C
libraries have to apply some common sense. That is a
fact (unfortunate perhaps) of the language.

And he did it because it was useful to write it
down somewhere. Don't quibble about where he wrote
it down.


> I think that you've painted yourself into a corner
> and you will be obliged to support this particular
> data representation and the implied algorithm forever.

Nobody is painting any corners, and nobody is
obliged to do anything. Chill out.


--
G. Jungman

Re: gsl Development Query

[E. Robert Tisdale]

Please take a look at
the Vector, Signal and Image Processing Library (VSIPL)

    http://www.vsipl.org/

Lots of really smart people have been working very hard on
a digital signal processing library API standard for over four years now.
The VSIPL ANSI C language binding is very similar
to the GSL ANSI C language binding.

Yes.  The GSL should allow application programmers
to create FFT objects like the "plan" object in FFTW.
The overhead of creating an FFT object is amortized
over several invocations of an FFT of a particular length.
The GSL really should include an entire "wing"
devoted to digital signal processing.  The application programmer
should be able to create and invoke convolver, correlation,
filter and decomposition objects as well as FFT objects.

Re: gsl Development Query

[Brian Gough]

E. Robert Tisdale writes:
 > The GSL should allow application programmers
 > to create FFT objects like the "plan" object in FFTW.
 > The overhead of creating an FFT object is amortized
 > over several invocations of an FFT of a particular length.

Yes, this is how it works -- see the documentation for
gsl_fft_complex_alloc in the manual for details.

PPC Linux Bugs, and a perplexing problem

[F J Franklin]

Discovered GSL only this week, and must say I'm *very* impressed! However...

> Until we have set up a separate bug reporting address, please report
> bugs to the GSL discussion list gsl-discuss@sourceware.cygnus.com

Compiled & checked happily enough on my (x86) desktop, but not on my (G3) 
laptop. (This is the first time I've ever submitted a bug report, so please 
forgive me if I have got the format all wrong.)

Has anyone succeeded in compiling & checking GSL-0.6 for PPC Linux?

Frank

Francis James Franklin
fjf@alinameridon.com

Diodorus the professor of logic died of shame because he could not at once 
solve a problem put to him in jest by Stilpo.
                                                       --- Pliny the Elder

=============================== Bugs ===============================

 (1) `make' Error
 (2) `make check' Errors
 (3) Puzzling PPC Problem

GSL 0.6

LinuxPPC 2000 on Apple iBook SE
kernel: linux-2.2.17pre13-ben1
compiler: gcc-2.95.2
libc & libm: glibc-2.1.3

(1) `make' Error
    ============

Neither of the following options is (or can be) set in config.h 
/* #undef HAVE_M68KLINUX_IEEE_INTERFACE */
/* #define HAVE_LINUX_IEEE_INTERFACE 1 */

FYI:

In /usr/include/fpu_control.h:
_FPU_SINGLE,_FPU_DOUBLE,_FPU_EXTENDED                     *Not* defined
_FPU_RC_NEAREST,_FPU_RC_DOWN,_FPU_RC_UP,_FPU_RC_ZERO      Defined
_FPU_MASK_IM,_FPU_MASK_ZM,_FPU_MASK_OM,_FPU_MASK_UM       Defined
_FPU_MASK_DM                                              *Not* defined
_FPU_MASK_PM                                              *Not* defined

However, defined in /usr/include/fpu_control.h is (cf. _FPU_MASK_PM ?)
/* masking of interrupts */
#define _FPU_MASK_XM  0x08 /* inexact */

(2) `make check' Errors
    ===================

(a) BLAS
    ====

Segmentation fault!
Source: In `int test_gbmv(void)' in `gsl-0.6/blas/test_blas_raw.c'

  gsl_blas_raw_dcopy(4, vector_4_d, 1, tmp_d, 1);
  gsl_blas_raw_dgbmv(CblasNoTrans, 4, 4, 1, 2, 2.0, matrix_gen_4_d, 4, vector_4_d, 1, 3.0, tmp_d, 1);
  s = ( tmp_d[0] != -10.0 || tmp_d[1] != -5.0 || tmp_d[2] != 7.0 || tmp_d[3] != 9.0 );
  gsl_test(s, "gsl_blas_raw_dgbmv A");
  status += s;

  gsl_blas_raw_dcopy(4, vector_4_d, 1, tmp_d, 1);
  gsl_blas_raw_dgbmv(CblasNoTrans, 4, 4, 1, 2, 2.0, matrix_gen_4_d, 4, vector_4_z, 2, 3.0, tmp_d, 1);
  s = ( tmp_d[0] != -10.0 || tmp_d[1] != -5.0 || tmp_d[2] != 7.0 || tmp_d[3] != 9.0 );
  gsl_test(s, "gsl_blas_raw_dgbmv B");
  status += s;

Commenting these out `removes' the problem... (more on this later; see below)

(b) MONTE [This is listed in KNOWN-PROBLEMS, I see...]
    =====

Testing double gaussian
FAIL: vegas(f2), dim=9, err=0.0003, chisq=0.6418 (0.499681001846611517 observed vs 1 expected)
FAIL: vegas_test

*** Rather suspicious factor of 2 here! (to state the obvious)

(c) SPECFUNC
    ========

FAIL:   gsl_sf_hyperg_1F1_int_impl(100, 400, 100.0, &r)
  test_hyperg.c 73
  expected:     907123037681.855
  obtained:    907123037681.8558    0.04068719885713332  4.4853e-14
  fracdiff: 4.709902361666654e-16
  value not within tolerance of expected value
     907123037681.855835   0.0406871988571333165
FAIL:   gsl_sf_hyperg_U_int_impl(-90, 1, 10, &r)
  test_hyperg.c 323
  expected: -1.898767714922189e+139
  obtained: -1.898767714922192e+139   1.543096725201747e+126  8.12683e-14
  fracdiff: 6.962999882451004e-16
  value not within tolerance of expected value
  -1.8987677149221916e+139  1.54309672520174719e+126
FAIL:   gsl_sf_hyperg_U_int_impl(-90, 10, 10, &r)
  test_hyperg.c 324
  expected: -5.682999988842067e+143
  obtained: -5.682999988842072e+143   4.618478923559706e+130  8.12683e-14
  fracdiff: 4.573953621928294e-16
  value not within tolerance of expected value
  -5.68299998884207181e+143  4.61847892355970615e+130
FAIL:   gsl_sf_hyperg_U_impl(1, 1.2, 2.0, &r)
  test_hyperg.c 360
  expected:   0.3835044780075603
  obtained:   0.3835044780075599   9.696010384367307e-16  2.52827e-15
  fracdiff: 5.066147605858526e-16
  value not within tolerance of expected value
    0.383504478007559879  9.69601038436730741e-16
FAIL:   gsl_sf_hyperg_U_impl(-50.5, 100.1, 40, &r)
  test_hyperg.c 387
  expected: 5.937463786613894e+91
  obtained: 5.937463786613886e+91   8.703956648586891e+79  1.46594e-12
  fracdiff: 6.703793421875711e-16
  value not within tolerance of expected value
  5.93746378661388571e+91  8.70395664858689101e+79
FAIL:   gsl_sf_hyperg_2F1_conj_impl(25, 25, 1, -0.5, &r)
  test_hyperg.c 424
  expected: 5.16969440956632e-06
  obtained: 5.169694409691514e-06   5.713800742253134e-17  1.10525e-11
  fracdiff: 1.210837176728465e-11
  value/expected not consistent within reported error
  5.16969440969151365e-06  5.71380074225313361e-17
FAIL:   gsl_sf_hyperg_2F1_conj_renorm_impl(9, 9, -1.5, -0.99, &r)
  test_hyperg.c 454
  expected:   0.1083402022947613
  obtained:   0.1083402022947615   4.247728659543883e-13  3.92073e-12
  fracdiff: 1.152850812762138e-15
  value not within tolerance of expected value
    0.108340202294761503  4.24772865954388333e-13
FAIL: Hypergeometric Functions

FAIL:   gsl_sf_laguerre_n_impl(90, 2.0, 100.0, &r)
  test_sf.c 873
  expected: -2.092904225954693e+20
  obtained: -2.09290422595469e+20      50333025.30505015  2.40494e-13
  fracdiff: 7.82835630833816e-16
  value not within tolerance of expected value
  -2.09290422595468952e+20     50333025.3050501496
FAIL: Laguerre Polynomials

FAIL:   gsl_sf_synchrotron_2_impl(0.01, &r)
  test_sf.c 1083
3
  obtained:   0.1083402022947615   4.247728659543883e-13  3.92073e-12
  fracdiff: 1.152850812762138e-15
  value not within tolerance of expected value
  expected:   0.2309807734222628
  obtained:   0.2309807734222623   1.877351636458429e-16  8.12774e-16
  fracdiff: 9.613120678191913e-16
  value/expected not consistent within reported error
    0.230980773422262337  1.87735163645842852e-16
FAIL: Synchrotron Functions

(3) Puzzling PPC Problem
    ====================

GSL-0.6 compiles and checks for x86-Linux, but not for iBook Linux. As far as 
I can tell, the code itself is not at fault, but I have isolated the problem 
in the following extract-example program. The output is given below; note how 
args 12 & 13 get messed up.

Can anybody tell me what, if anything, may be wrong with the following code?
(The program behaves as expected when compiled on my x86, so the problem is
probably PPC-specific.)

Interestingly, the problem vanishes if every instance of `double' is replaced
by `float'...

=============================== test.c ===============================

#include <stdio.h>

typedef enum CBLAS_TRANSPOSE
{  CblasNoTrans   = 111,
   CblasTrans     = 112,
   CblasConjTrans = 113
} CBLAS_TRANSPOSE_t;

const double     vector_4_d[] = { -2.0, -1.0, 0.0, 3.0 };
const double matrix_gen_4_d[] = 
{
  1.0,  0.0, -2.0, 1.0, 
  0.5,  3.0,  5.0, 1.0,
  2.0, -2.0, -1.0, 0.5,
  1.0, -1.0,  0.5, 0.0
};

void dgbmv (CBLAS_TRANSPOSE_t TransA,
            size_t M, size_t N, size_t KL, size_t KU,
            double alpha,
            const double A[], int lda,
            const double X[], size_t incX,
            double beta,
            double Y[], size_t incY)
{
  fprintf (stderr,"dgbmv: TransA = %d\n", (int) TransA);
  fprintf (stderr,"dgbmv: M      = %d\n", (int) M);
  fprintf (stderr,"dgbmv: N      = %d\n", (int) N);
  fprintf (stderr,"dgbmv: KL     = %d\n", (int) KL);
  fprintf (stderr,"dgbmv: KU     = %d\n", (int) KU);
  fprintf (stderr,"dgbmv: alpha  = %lf\n",alpha);
  fprintf (stderr,"dgbmv: A      = %p\n", A);
  fprintf (stderr,"dgbmv: lda    = %d\n", lda);
  fprintf (stderr,"dgbmv: X      = %p\n", X);
  fprintf (stderr,"dgbmv: incX   = %d\n", (int) incX);
  fprintf (stderr,"dgbmv: beta   = %lf\n",beta);
  fprintf (stderr,"dgbmv: Y      = %p\n", Y);
  fprintf (stderr,"dgbmv: incY   = %d\n", (int) incY);
}

int main()
{
  double tmp_d[32];

  fprintf (stderr,"dgbmv: tmp_d          = %p\n",tmp_d);
  fprintf (stderr,"dgbmv: vector_4_d     = %p\n",vector_4_d);
  fprintf (stderr,"dgbmv: matrix_gen_4_d = %p\n",matrix_gen_4_d);

  dgbmv(CblasNoTrans,4,4,1,2,2.0,matrix_gen_4_d,4,vector_4_d,1,3.0,tmp_d,1);
}

=============================== PPC output ===============================

dgbmv: tmp_d          = 0x7ffffb18
dgbmv: vector_4_d     = 0x10000730
dgbmv: matrix_gen_4_d = 0x10000750
dgbmv: TransA = 111
dgbmv: M      = 4
dgbmv: N      = 4
dgbmv: KL     = 1
dgbmv: KU     = 2
dgbmv: alpha  = 2.000000
dgbmv: A      = 0x10000750
dgbmv: lda    = 4
dgbmv: X      = 0x10000730
dgbmv: incX   = 1
dgbmv: beta   = 3.000000
dgbmv: Y      = 0x1
dgbmv: incY   = 805388944

=============================== x86 output ===============================

dgbmv: tmp_d          = 0xbffffad8
dgbmv: vector_4_d     = 0x8048668
dgbmv: matrix_gen_4_d = 0x8048688
dgbmv: TransA = 111
dgbmv: M      = 4
dgbmv: N      = 4
dgbmv: KL     = 1
dgbmv: KU     = 2
dgbmv: alpha  = 2.000000
dgbmv: A      = 0x8048688
dgbmv: lda    = 4
dgbmv: X      = 0x8048668
dgbmv: incX   = 1
dgbmv: beta   = 3.000000
dgbmv: Y      = 0xbffffad8
dgbmv: incY   = 1

gsl Development Query

[F J Franklin]

Dear All,

Is anybody working on multidimensional fourier transforms? I wouldn't mind
trying my hand at it, having had some (as in `a little') experience in
this area.

Q: Is there some kind of guide for developers?

I am particularly interested in issues of function-naming,
multi-dimensional array representation, inline vs workspace algorithms,
`ethics' of dynamic memory allocation, IEEE issues & testing/checking,...

Any and all comments appreciated...

Frank

Francis James Franklin
fjf@alinameridon.com

Diodorus the professor of logic died of shame because he could not at once 
solve a problem put to him in jest by Stilpo.
                                                       --- Pliny the Elder

Re: gsl Development Query

[Thomas Walter]

>>>>> "F" == F J Franklin <MEP95JFF@sheffield.ac.uk> writes:

    F> Dear All,
    F> Is anybody working on multidimensional fourier transforms? I wouldn't mind
    F> trying my hand at it, having had some (as in `a little') experience in
    F> this area.

    F> Q: Is there some kind of guide for developers?

    F> I am particularly interested in issues of function-naming,
    F> multi-dimensional array representation, inline vs workspace algorithms,
    F> `ethics' of dynamic memory allocation, IEEE issues & testing/checking,...

    F> Any and all comments appreciated...

    F> Frank

Hello, have a look at
    http://www.fftw.org

This is a fast fourier package for any size and multi dimensions.

Bye
Thomas


-- 
Was gibt sieben mal sieben?  Ganz feinen Sand. 8-)

----------------------------------------------
Dipl. Phys. Thomas Walter
Inst. f. Physiklische Chemie II
Egerlandstr. 3				Tel.: ++9131-85 27326 / 27330
91058 Erlangen, Germany			email: walter@pctc.chemie.uni-erlangen.de

Re: gsl Development Query

[Brian Gough]

F J Franklin writes:
 > Is anybody working on multidimensional fourier transforms? I wouldn't mind
 > trying my hand at it, having had some (as in `a little') experience in
 > this area.

Nobody is working on that currently, it is in the fft/TODO file
("simple multidimensional fft") and would be nice to have.  For
serious fft work we recommend using FFTW of course.

 > Q: Is there some kind of guide for developers?

There are some notes in CVS, doc/gsl-design.texi.  See
http://sources.redhat.com/gsl/ for details on anonymous CVS access.

Re: gsl Development Query

[Gerard Jungman]

Brian Gough wrote:
>
> For serious fft work we recommend using FFTW of course.

I agree. But isn't this an odd statement?
I really don't get it sometimes. If using FFTW
is the right thing to do, then why do we provide
our own implementation, with no explicit support for
connecting to outside implementations? This is the
worst of both worlds: tell people to use the other
thing, but don't provide any explicit way for them to
connect code using gsl abstractions to whatever
the external thing is.

And what exactly constitutes "serious work"?
Is GSL not for "serious work"?

There _should_ be gsl FFTs, for all dimensions,
but only in the sense that there should be
gsl defined _interfaces_ to such functionality.
These should be general enough to accomodate
any reasonable underlying implementation,
and should be a fixed target. It should be
possible for the user as well as the developer
to swap the underlying implementation at will.

Whether or not GSL itself provides a "native"
implementation is really a second order issue.
First design the interfaces, then connect them
to something that works (anything, to start),
and move on.

The same issue has arisen with the linear algebra
stuff. The world doesn't need another C implementation
of an FFT library any more than it needs another
implementation of BLAS and/or LAPACK, in any language.
I see now that the decision to provide a native GSL
implementation of linear algebra was probably a mistake.
I blame myself for that. But there is still time to
correct the problem there. Let's not make the same
mistake with FFTs.

If somebody wants to do something useful
(and we need all the help we can get), then
they should consider the following:

  1) We need a general definition of interfaces
     to FFT functionality, designed with GSL
     flavor, using GSL error reporting conventions,
     and probably paying at least some respect
     to the GSL data representations, conventions
     regarding stride, etc. This should be easy
     to do. All you have to do is look at available
     implementations, abstract out their common
     properties, and stir in a little GSL.

  2) We need a paradigm for connecting these GSL
     interfaces to external implementations. This
     must include the possibility of connecting
     to fortran implementations. Obviously, you
     should be able to connect to a "native"
     implementation as well (such as Brian's 1d
     implementation), if only as a fallback.

  3) Pick a default "example" implementation and
     make it as easy as possible for people writing
     GSL code to swap in that implementation. I tried
     to do this with BLAS by defining a set of GSL
     interfaces (the "gsl_blas_raw" interfaces) for
     wrapping external blas implementations and then
     providing a native implementation expressing
     those interfaces. However, maybe there is
     something wrong with this, since one of the
     most common messages on the mailing list is
     "what are all these undefined blas symbols?"
     People don't understand that they have been
     given the freedom to swap in any conformant
     implementation at link time, and so they
     have to choose. Frankly, I think some of
     the developers still don't get it (sigh),
     although it doesn't seem like such a bad
     idea to me.


>  > Q: Is there some kind of guide for developers?

By the way, I have been told that one of the
statements in the design document is something
about "no mixed language development". Some of the
developers seem to think that this implies we
should not provide interfaces to underlying
fortran implementations, and that we should
certainly not package up a wrapped fortran
implementation for distribution as a default
backend.

This sort of thinking is manifestly confused.
Is one of the design goals of GSL to waste
development time reimplementing perfectly
adequate and tested functionality from other
libraries? A fortran subroutine can be a
perfecty adequate tool for multiplying
matrices (after all, that is what fortran
does best).

Anyway, this issue of "mixed language"
support applies to the question of
how to provide linear algebra. But in
the case of FFTs, given the existence
of a high quality FFT implementation
in C (FFTW), I hope the correct path
is clear to all.


-- 
G. Jungman

Re: gsl Development Query

[Brian Gough]

What you've outlined is something much bigger that GSL (i.e. it's not
GSL).  GSL is just a collection of routines for numerical computing,
written in C with modern coding conventions.  

Maybe a "universal interface/library" of the type you suggest is
worthwhile -- but it would be a different project from GSL.

If you need wrappers for fortran codes they are available in other
places (e.g. Octave's library, NumPy, ...).  The idea of GSL is that
wrapping is never completely satisfactory and so it's better to
rewrite things.  However much you wrap unreadable fortran, it's still
unreadable fortran.

Re: gsl Development Query

[Gerard Jungman]

Brian Gough wrote:
> 
> What you've outlined is something much bigger that GSL (i.e. it's not
> GSL).

No. What I have described is an efficient way to
obtain the required functionality in GSL. These
are implementation details. If you feel there
are some fundamental flaws in existing implementations,
then you are free to reimplement. I don't see the
flaws in existing (and Free!) implementations
such as BLAS and FFTW.

And how can it possibly matter to the end user? All they
want is something that works, and that works in the right way,
without making undesirable assumptions about usage patterns
or data representation. These important issues are _interface_ issues.
Their only relation to the elemental entity turning the crank
under the hood is to specify certain constraints that any
such implementation must satisfy. If you find a solution
to those constraints, then you have satisfied the design
criteria. If the solution happens to already exist, then
so much the better; you can go home early.


> GSL is just a collection of routines for numerical computing,
> written in C with modern coding conventions.

Yes, that's what it says in The Big Book of GSL. Now,
rather than quoting the precious scripture to me, why don't
you apply your native intelligence and think about whether or
not it makes any sense.

The so-called "design" of this library was never anything
more than some vague wish list. The original idea that you
could provide "loosely coupled" libraries for the different
tasks was wrong. This was made clear when Brian and I got
down to worrying about the various multi-dimensional
functionalities which depend explicitly on data structures
for vectors and matrices and linear algebra support. The
linear algebra bottleneck should be obvious to everyone
who has done any numerical computing; if you don't solve
this problem first, you are screwed.

When you take the time to look at the requirements for
implementing higher level algorithms, you realize
that the design has to be layered, and you have
to solve the fundamental problems first. These include
things like control of and access to the floating point
capabilities of the machine, data representation issues,
and basic linear algebra.

When you (Brian) and I faced these problems squarely, about
two years ago, we found it necessary to redesign and
reimplement large sections of the project, fixing
things related to the issues above, as well as
doing the right thing at higher levels, such as
introducing a consistent framework for iteration.

When you say "well, we're just a collection of tools",
you are just echoing the original design goals, which
have been demonstrated, by your own work, to be
confused.


> Maybe a "universal interface/library" of the type you suggest is
> worthwhile -- but it would be a different project from GSL.

I'm not talking about a "universal interface/library",
whatever that is. I'm talking about the right way to
provide an underlying implementation for some of the
required functionality.


> If you need wrappers for fortran codes they are available in other
> places (e.g. Octave's library, NumPy, ...).

Yes, they are. So?

Nobody "needs" wrappers for fortran codes. They need
tools that do certain tasks. Wrapping fortran just happens
to be a convenient way to get certain functionality.
That's why people do it, not because they like to
know that a little fortran monkey is turning the
crank in their greasy little simulated worlds.
Who cares what turns the crank, as long
as it turns?


> The idea of GSL is that
> wrapping is never completely satisfactory and so it's better to
> rewrite things.

Let me paraphrase this for the general reader: The idea of
GSL is to write half-assed implementations of some ill-defined
collection of numerical foo, sometimes of functionality
that already exists, with no real thought going into the
overall design, and to pat each other on the back for
having done everything the hard way and put it under
the GPL. Well, GPL'ed foo is still foo; that's one thing
I know for sure after having watched the free software
world work for the last ten years.

You have this idea that GSL is just some thing that
we are tinkering with in the backyard, not really
"important" enough to do the right way; after all,
nobody's going to the moon riding on a GSL rocket.
This is insulting to everyone involved, including
yourself.


> However much you wrap unreadable fortran, it's still
> unreadable fortran.

Huh? What does this mean? This is a non sequitur.


-- 
G. Jungman

Re: gsl Development Query

[Randall Judd]

Hi Group,

     Just a couple of comments from somebody who does not know very much
about gsl. I work on VSIPL and started lurking about to see how well gsl
and VSIPL will interface, since gsl is working on a lot of stuff VSIPL does
not do.

     Unlike gsl, VSIPL hides most stuff. We are a specification, not a
library. Implementations are free to implement blocks and views any way
they want as long as they meet the specification. Part of the specification
is the interface which allows people to bring data into or export data from
VSIPL.

     Now I have also been working on a VSIPL implementation, which is open
source, free and written in ANSI C, and available at the vsipl web site.
Long ago I needed an FFT and used FFTW in my implementation. Since the FFTW
copyright is too strong for the purposes of my implementation I just
encapsulated it and directed people to go to the FFTW site, get that
library, and then link it in when they built mine.

      The main problem with FFTW is 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). However it is not "comfortable" just assuming
that the two will work together. So eventually I wrote my own FFT stuff
(which is pretty ugly since I don't know much about FFT's) just to make the
FFT fit better with the VSIPL spec. 

      I guess my point is you probably should have your own FFT's. It may
not be reasonable to make an interface to FFTW and keep the other
interfaces in your library coherent. Frequently the FFT speed is not that
important and, as long as your FFT speed is reasonable, it may be fine for
most peoples use. I assume your fft inputs and outputs will fit well with
other functions you have defined and this may be more important to your
users. 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. 

                  Randy Judd

Re: gsl Development Query

[Gerard Jungman]

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?

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


> I guess my point is 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 blas. All you have to do is specify an appropriate middle
layer in each case.

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.


Thanks for your comments.

-- 
G. Jungman

Re: gsl Development Query

[Randall Judd]

Well, Gerard gave me a really long answer/comment and I appreciate it. I
may have some reply on other parts of it in a few days but for now I will
stick with the first couple of paragraphs since that is what interests me
mostly.

First I would like to point out some differences between what I perceive
you are doing and what VSIPL did. VSIPL is a specification. The intent from
the beginning was to get multiple vendors to supply their own optimized
VSIPL libraries for their own product. What I perceive GSL is doing is a
single library, not a specification. You don't expect anybody else to take
your document and write another version of GSL, at least I don't think you do.

So when VSIPL meets we have several vendors who participate and are
interested in not doing anything they don't need to in order to keep their
implementation cost down. So a lot of functionality you have in GSL was not
put in VSIPL. So I guess what I want is to have a clean interface between
VSIPL and GSL so that I can recover some of the functionality that the
vendors did not want. 

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

Just a warning. I don't really know very much about C++ and I like C a lot.
I can see a need to learn more about C++ and from time to time make an
effort, but so far I have not been convinced it is better than C. When I
actually learn it no doubt I will come around.

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

Well I am not sure what you are talking about below, so I will go into some
detail as to how  VSIPL would interface to GSL. Right now I think VSIPL
should interface fine to GSL, so hopefully if you understand a little more
you wont change something to screw it up.
 
>By the way, I looked at the VSIPL draft to see what
VSIPL 1.0 is actually final now. It is available on the Documents page of
the VSIPL site. On that page you will find something called, I think, "The
basics requirements document". You might want to read that instead. It is
short, and covers most of how VSIPL is supposed to work. Depending on when
you read the Draft VSIPL may have changed a lot from what you read.

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

First of all in VSIPL a view is always bound to a block. The blocks in GSL
just seem to be a place to keep your memory storage until your ready to
destroy it. In VSIPL blocks are more important. To find the begining of
your view in GSL you store a pointer to the beginning of data. In VSIPL we
store an offset into the block. The vendor may do other stuff under the
covers, but what the user sees is an offset. 

How a block in VSIPL is implemented is vendor dependent, and the blocks are
incomplete typedefs. For instance what a user sees in the vsipl header file
is something like

struct vsip_blockobject_f;
typedef struct vsip_blockobject_f vsip_block_f;

for a real block of type float. The vendor has a complete type definition,
of course, but he does not give that to the user. What the user gets is a
set of VSIPL defined functions to create and manipulate the block. "Views"
(vector views, matrix views) on the block are treated the same way.

struct vsip_vviewobject_f;
typedef struct vsip_vviewobject_f vsip_vview_f;

So to create a vector you could do

vsip_block_f *a_block = vsip_blockcreate_f(N,VSIP_MEM_NONE);
vsip_vview_f *a_view = vsip_vbind_f(a_block,0,1,N);

So a_view is a vector of length N, stride 1, and offset into the block of 0.

Since a_view is nice and compact we could have done the same thing with a
convenience function
vsip_vview_f *a_view = vsip_vcreate_f(N,VSIP_MEM_NONE);

Now the problem is that all this stuff is opaque. Vendors want to own the
data. Memory management was a big problem for the forum in the beginning,
and the fix was to abstract memory away from the hardware by using the
block concept along with abstract data types and incomplete typedefs.

So how does the user get data in and out of VSIPL? Well we defined some
special functions to do this. By the way, for the stuff above there is no
legal way in VSIPL for a user to get a pointer to the actual data storage
for the block.

What VSIPL does is define a blockbind function. What blockbind does is to
associate some regular memory with a VSIPL block. Then both VSIPL and the
user can get to that memory. However VSIPL insists on owning any memory it
might use, so we added functions to admit the block (and memory) to VSIPL
and a release function to give ownership of the data back to the user.
Admit and release are important. VSIPL does not necessarily use the user
memory in the expected way. It may not use it at all. The purpose of admit
and release is to force data consistency between the block and the memory
during a state transition.

So if you wanted a gsl vector, and a vsipl vector on the same data space
you could do (I don't know much about gsl so there may be an error here)

gsl_vector *a_gsl_vview = gsl_vector_alloc(N);
vsip_block_d *a_vsipl_block = vsip_blockbind_d(
        a_gsl_vview->block->data, N, VSIP_MEM_NONE);
vsip_vview_d *a_vsipl_vview = vsip_vbind_d(a_vsipl_block, 0, 1 , N);

To make a_vsipl_vview usable by VSIP (and make sure the data agrees with
a_gsl_vview) I would do

vsip_blockadmit_f(a_vsipl_block,VSIP_TRUE);

Once VSIPL is done and I want to do something with gsl (I should not use
the GSL vector while the associated data is admited) I would do
vsip_blockrelease_f(a_vsipl_block,VSIP_TRUE);


NOTE: It is not a mistake when I attach VSIPL to the data in the GSL block
instead of the GSL view. You don't want the same data segment attached to
two different VSIPL blocks. This is one of the problems I have with the way
GSL did this, but it appears that you can with a little care create
matching gsl and VSIPL views for any number of subviews on the original
data set. VSIPL uses an offset to do this. VSIPL views are more flexible
than GSL views, so not any VSIPL view is mapable to a GSL view, but with
knowledge of both libraries it seems to be pretty reasonable to use them
together, at least at first glance.

So, by now your probably really confused. I just hope I wasn't. If you have
any questions let me know.

                 Randy Judd.
       

Re: gsl Development Query

[Gerard Jungman]

Randall Judd wrote:
> 
> VSIPL 1.0 is actually final now. It is available on the Documents page of
> the VSIPL site. On that page you will find something called, I think, "The
> basics requirements document". You might want to read that instead.

Yes, I just read it. It's marked "DRAFT 0.6", but I suppose I
should ignore that. There are also the "Core" and "Core Lite" documents,
which seem to be just tables of types and function names.

Anyway, as you say, it certainly seems like there is
nothing to prevent you from interfacing with GSL in
a natural way. I'm not even sure what we could do
to screw that up. We're just a little closer to
the underlying memory, which, I think, is natural
for our purposes.


> First of all in VSIPL a view is always bound to a block. The blocks in GSL
> just seem to be a place to keep your memory storage until your ready to
> destroy it. In VSIPL blocks are more important. To find the begining of
> your view in GSL you store a pointer to the beginning of data. In VSIPL we
> store an offset into the block. The vendor may do other stuff under the
> covers, but what the user sees is an offset.

I understand. Your blocks are already at a higher level
of abstraction, whereas ours really serve at the lowest
level, for basic heap management. Again, no problems.
I think.


> The vendor has a complete type definition,
> of course, but he does not give that to the user. What the user gets is a
> set of VSIPL defined functions to create and manipulate the block. "Views"
> (vector views, matrix views) on the block are treated the same way.

Maybe I should say explicitly that I don't think we should
expose things like the gsl_block attributes directly either.
The only sense in which we do is that the struct is there
in the header file. But I have always felt that this is
just a matter of style in C; for instance, I can go look
up lots of things in the standard C header files which I would
never rely on in an application program.

As you say, we're writing a library, not a specification.
So the struct defs have to go somewhere. The only reason
I make this obvious statement is because it seems to be
one of those things that Tisdale is ranting about (at
least as far as I understand anything he says).


> What VSIPL does is define a blockbind function. What blockbind does is to
> associate some regular memory with a VSIPL block. Then both VSIPL and the
> user can get to that memory. However VSIPL insists on owning any memory it
> might use, so we added functions to admit the block (and memory) to VSIPL
> and a release function to give ownership of the data back to the user.
> Admit and release are important. VSIPL does not necessarily use the user
> memory in the expected way. It may not use it at all. The purpose of admit
> and release is to force data consistency between the block and the memory
> during a state transition.

Ok.


> So if you wanted a gsl vector, and a vsipl vector on the same data space
> you could do (I don't know much about gsl so there may be an error here)
> 
> gsl_vector *a_gsl_vview = gsl_vector_alloc(N);
> vsip_block_d *a_vsipl_block = vsip_blockbind_d(
>         a_gsl_vview->block->data, N, VSIP_MEM_NONE);
> vsip_vview_d *a_vsipl_vview = vsip_vbind_d(a_vsipl_block, 0, 1 , N);

A minor point, but maybe important. Why not use
the gsl_vector_ptr() method to get the pointer
to the underlying type, rather than "a_gsl_vview->block->data"?
Anyway, that is supposed to provide the necessary abstraction
away from the explicit struct attributes.
[Ahah. I see that this is the point, from your
 comments below. Ok, let's move on to that then...]

[Aside to Brian: Is the damn definition for gsl_vector_ptr() missing?
 I can't find the implementation anywhere...]


> NOTE: It is not a mistake when I attach VSIPL to the data in the GSL block
> instead of the GSL view. You don't want the same data segment attached to
> two different VSIPL blocks. This is one of the problems I have with the way
> GSL did this, but it appears that you can with a little care create
> matching gsl and VSIPL views for any number of subviews on the original
> data set.

Ok. This is the part I need to understand. Let me state what I
understand and you can tell me if I am missing the point. The problem
is that you would like to attach VSIPL blocks to GSL views,
but the semantics do not match because GSL views might overlap,
and VSIPL blocks cannot. So instead you have to go directly to
the GSL block, which breaks any encapsulation we might have hoped for.

If this is the problem, then the right solution is, as you do,
to go directly to the block. So then my question is, what
would you like to see us do with gsl_block in order to make
this mapping easier?

Do you think gsl needs another abstraction in the middle,
which is closer to the VSIPL block notion? That might be
overkill, but it wouldn't be too hard.

Is the real question here ownership and the possibility
of overlaps, or is it something else?


> VSIPL uses an offset to do this. VSIPL views are more flexible
> than GSL views, so not any VSIPL view is mapable to a GSL view, but with
> knowledge of both libraries it seems to be pretty reasonable to use them
> together, at least at first glance.

GSL views may be less flexible simply because we have not
fleshed them out more. But I don't see them as fundamentally
less flexible. Do you agree, or do you think there is a
real flaw there?

Anyway, if I understand what you mean, then the right thing
to do is to make it easier to do the correct
mapping, which is blocks to blocks, with the 
two separate view abstractions remaining in
their own worlds. There may be no reason to
attempt to map the two different views to
each other in any way.


> So, by now your probably really confused.

Probably. You can judge.


Thanks.

-- 
G. Jungman

Re: gsl Development Query

[Randall Judd]

At 03:58 PM 09/02/2000 -0600, Gerard Jungman wrote:
>Randall Judd wrote:
>> 
>> VSIPL 1.0 is actually final now. It is available on the Documents page of
>> the VSIPL site. On that page you will find something called, I think, "The
>> basics requirements document". You might want to read that instead.
>
>Yes, I just read it. It's marked "DRAFT 0.6", but I suppose I
>should ignore that. There are also the "Core" and "Core Lite" documents,
>which seem to be just tables of types and function names.
I will check it out and complain to the proper person. I wrote it
originally, but then I lost control.

The Core, and Core Lite documents are what are called profiles. No vendor,
except perhaps SKY, has completed all of VSIPL yet. But they wanted to be
able to say they are VSIPL compliant. So we defined profiles to give them a
starting point and a way to claim compliance to something. Core Lite is
only about 129 function and is a very easy starting point for most vendors.
A great many, surprisingly so, current Navy applications appear to be
mappable to just a core lite implementation. The core profile has about 500
or so functions and include more complicated things. In the future their
may be more profiles. A profile is basically a defined set of functions
that has been approved by the VSIPL forum.

>
>Anyway, as you say, it certainly seems like there is
>nothing to prevent you from interfacing with GSL in
>a natural way. I'm not even sure what we could do
>to screw that up. We're just a little closer to
>the underlying memory, which, I think, is natural
>for our purposes.
I think I agree with this. I don't necessarily think you did it the best
way though. I don't know what all design constraints the GSL forum placed
on themselves.

>
>
>> First of all in VSIPL a view is always bound to a block. The blocks in GSL
>> just seem to be a place to keep your memory storage until your ready to
>> destroy it. In VSIPL blocks are more important. To find the begining of
>> your view in GSL you store a pointer to the beginning of data. In VSIPL we
>> store an offset into the block. The vendor may do other stuff under the
>> covers, but what the user sees is an offset.
>
>I understand. Your blocks are already at a higher level
>of abstraction, whereas ours really serve at the lowest
>level, for basic heap management. Again, no problems.
>I think.
>
>
>> The vendor has a complete type definition,
>> of course, but he does not give that to the user. What the user gets is a
>> set of VSIPL defined functions to create and manipulate the block. "Views"
>> (vector views, matrix views) on the block are treated the same way.
>
>Maybe I should say explicitly that I don't think we should
>expose things like the gsl_block attributes directly either.
>The only sense in which we do is that the struct is there
>in the header file. But I have always felt that this is
>just a matter of style in C; for instance, I can go look
>up lots of things in the standard C header files which I would
>never rely on in an application program.
>
>As you say, we're writing a library, not a specification.
>So the struct defs have to go somewhere. The only reason
>I make this obvious statement is because it seems to be
>one of those things that Tisdale is ranting about (at
>least as far as I understand anything he says).
>
>
>> What VSIPL does is define a blockbind function. What blockbind does is to
>> associate some regular memory with a VSIPL block. Then both VSIPL and the
>> user can get to that memory. However VSIPL insists on owning any memory it
>> might use, so we added functions to admit the block (and memory) to VSIPL
>> and a release function to give ownership of the data back to the user.
>> Admit and release are important. VSIPL does not necessarily use the user
>> memory in the expected way. It may not use it at all. The purpose of admit
>> and release is to force data consistency between the block and the memory
>> during a state transition.
>
>Ok.
>
>
>> So if you wanted a gsl vector, and a vsipl vector on the same data space
>> you could do (I don't know much about gsl so there may be an error here)
>> 
>> gsl_vector *a_gsl_vview = gsl_vector_alloc(N);
>> vsip_block_d *a_vsipl_block = vsip_blockbind_d(
>>         a_gsl_vview->block->data, N, VSIP_MEM_NONE);
>> vsip_vview_d *a_vsipl_vview = vsip_vbind_d(a_vsipl_block, 0, 1 , N);
>
>A minor point, but maybe important. Why not use
>the gsl_vector_ptr() method to get the pointer
>to the underlying type, rather than "a_gsl_vview->block->data"?
>Anyway, that is supposed to provide the necessary abstraction
>away from the explicit struct attributes.
I guess I don't understand why this is necessary?

>[Ahah. I see that this is the point, from your
> comments below. Ok, let's move on to that then...]
>
>[Aside to Brian: Is the damn definition for gsl_vector_ptr() missing?
> I can't find the implementation anywhere...]
>
>
>> NOTE: It is not a mistake when I attach VSIPL to the data in the GSL block
>> instead of the GSL view. You don't want the same data segment attached to
>> two different VSIPL blocks. This is one of the problems I have with the way
>> GSL did this, but it appears that you can with a little care create
>> matching gsl and VSIPL views for any number of subviews on the original
>> data set.
>
>Ok. This is the part I need to understand. Let me state what I
>understand and you can tell me if I am missing the point. The problem
>is that you would like to attach VSIPL blocks to GSL views,
>but the semantics do not match because GSL views might overlap,
>and VSIPL blocks cannot. So instead you have to go directly to
>the GSL block, which breaks any encapsulation we might have hoped for.
>

Ok. Here is the main problem in your understanding. I don't attach VSIPL
blocks to GSL views. The idea is to associate (a word I prefer to attach)
VSIPL blocks with GSL blocks (not views). Then define VSIPL vector and/or
matrix views which map the block data space the same way that the GSL views
map the GSL blocks data space. 

As I understand it your GSL block contains some data space. Then you slice
this data space up. The first of your vector views creates the block, data
space, and a vector. The vector view contains a pointer to the block, plus
a pointer to the beginning of the data. When you take of subview of the
vector you set the block pointer to NULL, get a pointer to the beginning of
the subviews data space, and then set the attributes of the subview to the
proper stride and length for the new vector, which must not exceed the data
space of the block.

In VSIPL every view contains a pointer to the block. The block contains the
entire possible data space. For a GSL vector view attached to a block the
view contains a pointer to the beginning of data, a length equal to the
data space length, and a stride equal to one. The corresponding VSIPL
vector view contains an offset of zero, a stride of one, and a length equal
to the length of the block.

So when I create a block I would want the block to encompass the entire GSL
block.

Now I could create a new block to do the GSL subview, but the new block
would encompass some of the same data space as the previous VSIPL block.
This would be an error from VSIPL's point of view. What I do in VSIPL is
set an offset into the data space so my subview points to the same data
space as the GSL subview. Then the stride and length would be the same as GSL.

For example

 gsl_vector *a_gsl_vview = gsl_vector_alloc(N);
/* we assume N >> 4 */
 gsl_vector *a_gsl_subview = gsl_vector_subvector(a_gsl_vview, 4, N/2);
 gsl_vector *another_gsl_subview = gsl_vector_subviector(a_gsl_subview,2,N/4);
 vsip_block_d *a_vsipl_block = vsip_blockbind_d(
         a_gsl_vview->block->data, N, VSIP_MEM_NONE);
 vsip_vview_d *a_vsipl_vview = vsip_vbind_d(a_vsipl_block, 0, 1 , N);
 vsip_vview_d *a_vsipl_subview = vsip_vbind_d(a_vsipl_block,4,1, N/2);
 vsip_vview_d *another_vsipl_subview = vsip_vbind_d(a_vsipl_block,6,1, N/4);

/* or */
/* vsip_vview_d *a_vsipl_subview = vsip_vsubview_d(a_vsipl_vview,4,N/2); */
/* vsip_vview_d *another_vsipl_subview = vsip_vsubview_d(a
vsipl_sub_view,2,N/4); */

>If this is the problem, then the right solution is, as you do,
>to go directly to the block. So then my question is, what
>would you like to see us do with gsl_block in order to make
>this mapping easier?
Well your GSL block is fine, as far as I can tell. I don't know how you
have done complex blocks. The problem is the view. It is not a big problem,
unless for some reason you are trying to back the information out instead
of knowing it up front. If you supplied an offsets instead of a data
pointer I would have an easy map between VSIPL and GSL. Right now I have to
know the offset. I can't get it from the view. If I know the parent view I
can figure it out, but it is not something one can get from the view
itself. If you kept the block around it would be doable, but you seem to
keep the block in only one view. I think you do this to control when the
data is destroyed, but their are better ways to do this. I can't think of
any STRONG reason you need to change anything. It is not what I would have
done, but it works. 

Actually, if you just added an offset and length to the view, so the
current pointer minus the offset would point to the beginning of the block,
and the length would be the length of the original block, then you would
have all the information of the original block carried with the new view
and it would not change any of your current programs. It might be a pain
keeping the offset information current, but a minor pain.

If I were doing it from scratch I probably would have made a set of very
tightly coupled blocks and views as abstract data types with functions to
return memory pointers, lengths and strides, and also to set memory
pointers lengths and strides. The internals of the blocks and views are not
important as long as your programers can get to the underlying memory
information for whatever particular function they are working on. This
would give you a support library which could be controlled by the main
forum, and all any of the other GSL function programers would need to do is
use this well understood support library when building their function. Of
course the initial support library definition may have been hard to define.
I know how these meetings go. But once done it would give you a common API
while allowing a diverse set of programmers to do their own thing. Right
now your API is not very consistent. 

>
>Do you think gsl needs another abstraction in the middle,
>which is closer to the VSIPL block notion? That might be
>overkill, but it wouldn't be too hard.
Since you are doing a C library, I don't see any reason to do this
abstraction. You may want to add methods to allow people to attach special
memory to one of you blocks, but since your block is already public they
can do that as is.

>
>Is the real question here ownership and the possibility
>of overlaps, or is it something else?
VSIPL did the memory abstraction of blocks for portability. Their is no way
to write portable software on diverse hardware if the user manages memory
directly. The VSIPL library manages memory, and a "block" looks the same on
any system, independent of hardware. One of the major drivers of VSIPL is a
need to port code from old hardware to new hardware. The cost of the code
port is huge. With the advent of COTS (commercial of the shelf) in
government programs when the upgrade cycle is a few years instead of
decades, the software port needs to be easy. The need for portable legacy
code is the driving factor behind VSIPL.

VSIPL has some of the same type of overlap requirements as GSL except we
allow a fair amount of in-place operations.

>
>
>> VSIPL uses an offset to do this. VSIPL views are more flexible
>> than GSL views, so not any VSIPL view is mapable to a GSL view, but with
>> knowledge of both libraries it seems to be pretty reasonable to use them
>> together, at least at first glance.
>
>GSL views may be less flexible simply because we have not
>fleshed them out more. But I don't see them as fundamentally
>less flexible. Do you agree, or do you think there is a
>real flaw there?
>
You define all your matrix views as being row major, unit stride in the row
direction. I assume when your people do their programing they use this
assumption to produce higher performance code, so I would consider it to be
fundamentally less flexible. Fleshing them out more won't fix all the
underlying code.

I don't know that this is a bad thing. It's just the way you did it.

VSIPL matrix views are either row major or column major, and we allow non
unit strides in the major direction. Also all views can be coerced into
looking at any part of the block. Their is no fundamental difference
between a parent view and a subview in VSIPL. You can make a subview a
clone of the parent by setting it's attributes to be the same as the parent. 

Also, I am not sure how you do complex exactly. I see a complex block, but
I have not found where you define it, or how a view is mapped on it. VSIPL
can create (real) views of the real or imaginary parts of a complex view
(this is a little tricky because a real view is bound to type vsip_block_f
and c complex view is bound to type vsip_cblock_f). You actually view the
same data space, and if you change the data in one view the corresponding
data will change in the other view.

I don't really know enough about how your views operate. For instance I see
where you can create a block, but I don't see what you can do with it once
it's created. It seems like, the way you have designed things, you would
always create a vector or matrix first. I suppose you can always brute
force the block into a vector. You have access to all your structures.
VSIPL has functions to handle all this.

>Anyway, if I understand what you mean, then the right thing
>to do is to make it easier to do the correct
>mapping, which is blocks to blocks, with the 
>two separate view abstractions remaining in
>their own worlds. There may be no reason to
>attempt to map the two different views to
>each other in any way.
This is right. Most people will probably not need to have a one to one
mapping of views between VSIPL and GSL. In fact they may want entirely
different views of the same data space, one in VSIPL and one in GSL. My
main goal here is to see how the two libraries can be used together to
allow users the most power to produce applications. What we don't know
won't help our efforts.

           Randy Judd

Re: gsl Development Query

[Brian Gough]

Randall Judd writes:
 > In VSIPL every view contains a pointer to the block. 

I think we can accomodate this in GSL by not setting the block field
to NULL for views, and having a separate field in the struct to
indicate ownership of the memory instead.  Then one would always be
able to determine the underlying block for any view.

This seems straightforward to support, so I will do that -- it only
needs a change in the struct defn and a few associated routines.

Re: gsl Development Query

[Randall Judd]

At 06:45 PM 09/04/2000 +0100, Brian Gough wrote:
>Randall Judd writes:
> > In VSIPL every view contains a pointer to the block. 
>
>I think we can accomodate this in GSL by not setting the block field
>to NULL for views, and having a separate field in the struct to
>indicate ownership of the memory instead.  Then one would always be
>able to determine the underlying block for any view.
>
>This seems straightforward to support, so I will do that -- it only
>needs a change in the struct defn and a few associated routines.
>
This would be good. I would then have enough information from the view to
create a matching VSIPL view and to check and see if the gsl block was
owned by the gsl view or not.

In VSIPL we leave it up to the user to keep track of what belongs with what
block and it is the users responsibility not to destroy the block before
everything is done using it. We have destroy (what your spec calls free)
functions for both views and blocks and an all destroy if we want to
destroy both view and block attached to view. For instance 
vsip_blockdestroy_f(vsip_block_f* vsip_vdestroy_f(vsip_vview_f*)) 
or 
vsip_valldestroy_f(vsip_vview_f*). 

For a developement mode VSIPL library we use a block attribute which
contains a number indicating valid when the block is made. When the block
is destroyed we set the number to invalid. This allows an implementation
(the user could not do this) to check and see if the block is valid or not.
Of course you don't check for invalid. The implementation needs to check
this in VSIPL because of the way we designed the library.

In section 6.1 of the GSL manual I downloaded I notice typdefs for complex
blocks, but I don't see anywhere where they are defined. Are the in the
manual? How is the stride handled? VSIPL strides are in terms of block
elements so the stride for real and the stride for complex are handled the
same. Historically many vector libraries do strides in terms of real
elements, even for complex. How does GSL handle this? Are there any
functions yet that need these comples blocks?

I run across an error from time to time in the manual. Is their a place to
send these too? For instance the example on page 246 has in the first for
loop,
REAL(z,i) = 0.0; which should be I think REAL(data,i)=0.0;

I spent (and am spending) a fair amount of time editing VSIPL specs and I
know how much trouble it is to get rid of errors. I just don't know if you
want me sending them to this mailing list.

               Randy Judd

Re: gsl Development Query

[Brian Gough]

Gerard Jungman writes:
 > No. What I have described is an efficient way to
 > obtain the required functionality in GSL.

My basic point: we don't need to put the functionality of FFTW in GSL.
People should just use FFTW if that is what they need. 

BLAS is a separate issue. If we didn't provide a sample BLAS
implementation then GSL would not install out-of-the-box, since BLAS
is used internally in the library. That's the only reason we provide
one.