PPC Linux Bugs, and a perplexing problem

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

Re: PPC Linux Bugs, and a perplexing problem

[Brian Gough]

Thanks for the bug report. 

1) Regarding the make error, ieee support for ppclinux has been added in
the CVS version of GSL.  See http://sources.redhat.com/gsl/ for
information on accessing it. 

2) Thanks for the listing of the make check errors -- we'll work on
those for the next release.

3) Regarding BLAS and the corrupted arguments -- hopefully someone can
explain that, it looks like a compiler problem. If so, a simple
function like void foo (int a, int b, int c, ...) might show the same
behavior for a sufficiently large number of arguments.

best regards

Brian Gough

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

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

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: PPC Linux Bugs, and a perplexing problem

[F J Franklin]

On Sat, 26 Aug 2000, Brian Gough wrote:

> 3) Regarding BLAS and the corrupted arguments -- hopefully someone can
> explain that, it looks like a compiler problem. If so, a simple
> function like void foo (int a, int b, int c, ...) might show the same
> behavior for a sufficiently large number of arguments.

Okay, turns out to be a gcc compiler bug, fixed in later versions of gcc.
If anyone's interested, the following code tests for it...

Frank

=============================================================================
dnl Early versions of gcc-2.95.2 for PPC have a bug when passing lots
dnl of arguments; test for this...
dnl (This test assumes cross-compiling is okay: ...? )

if test x"$ac_cv_prog_gcc" = xyes; then
 AC_MSG_CHECKING([for PPC gcc bug])

 AC_TRY_RUN([
  #include<stdlib.h>
  int test4bug (void* v1,void* v2,void* v3,void* v4,void* v5,void* v6,void* v7,
                double* ptr_1,void* v8,double d,double* ptr_2)
  { return ((ptr_1 == ptr_2) ? 0 : 1); }
  int main ()
  { double ptr[] = {0}; exit (test4bug (0,0,0,0,0,0,0,ptr,0,0.0,ptr)); }
 ],
 ac_ppc_gcc_bug=no,
 ac_ppc_gcc_bug=yes,
 ac_ppc_gcc_bug=no)

 if test $ac_ppc_gcc_bug = yes; then
  AC_MSG_RESULT(yes)
  AC_MSG_ERROR([gsl will not compile correctly; please update your compiler])
 else
  AC_MSG_RESULT(no)
 fi
fi
=============================================================================
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: PPC Linux Bugs, and a perplexing problem

[Brian Gough]

Thanks. It's good that you found the cause of the problem. I'll will
add the test and note it in the KNOWN-PROBLEMS file to save trouble
for other ppc users.