GSL

GSL

[Jens Lenge]

Hi Redhat team,

I heard about the GSL project and have some questions about it.
I hope you can help me.

GSL is often said to become a complete rewrite of LAPACK in native C, thus
offering LAPACK functionality without depending on F2C or any underlying
FORTRAN libs (like CLAPACK does).

==> Is it the goal of GSL to build a "native C" equivalent of LAPACK
functionality?
==> Will it depend on (or include) an underlying BLAS library?
==> If so, is that BLAS also native C?
==> Who is actually doing GSL? Is Redhat involved?
==> When used as a basis for a C++ wrapper lib, will the GSL data structures
be compatible with STL containers or the C++ complex<double> class?

Regards, Jens

--
Dipl.-Ing. Jens Lenge
Universitaet Dortmund, ET-HFT
lenge@hft.e-technik.uni-dortmund.de

Re: GSL

[Gerard Jungman]

Hi.

Jens Lenge wrote:
>
> GSL is often said to become a complete rewrite of LAPACK in native C, thus
> offering LAPACK functionality without depending on F2C or any underlying
> FORTRAN libs (like CLAPACK does).

Hmmm...


> ==> Is it the goal of GSL to build a "native C" equivalent of LAPACK
> functionality?

No. That would be too much. And we really don't see the
point, for several reasons. First, if what you want is
a C interface to LAPACK, then CLAPACK is probably more
appropriate. Second, effort on new native implementations
of this sort of thing would be better spent on C++, where
a lot of recent work has made great progress.

GSL has a native vector/matrix implementation, and it
would be really good to have a LAPACK interface for
working with these containers. But it could be only
a wrapper over CLAPACK. The amount of work required
in reimplementing LAPACK makes the idea insane.


> ==> Will it depend on (or include) an underlying BLAS library?

In my view, it should depend on a layer over a BLAS implementation.
The feature of the GSL work is that it is really just a layer
over a BLAS implementation, so in principle you can plug in
any implementation you like. The fact that GSL also comes
with a (incomplete) BLAS implementation of its own, I consider
to be of secondary importance.


> ==> If so, is that BLAS also native C?

As mentioned above, there is an incomplete native C BLAS
included in GSL. I wrote it. I think it was a bad idea,
but I was swayed by all the talk about having a "clean",
"new", "native C" implementation. Well, it's "new" anyway,
and it is written in C ...


> ==> Who is actually doing GSL? Is Redhat involved?

There is a link to the AUTHORS file on the web page, where
you will find all those people who have contributed substantial
code to the project. However, it's fair to say that the bulk of
the work has been done by Brian Gough and me.


> ==> When used as a basis for a C++ wrapper lib, will the GSL data structures
> be compatible with STL containers or the C++ complex<double> class?

Well, that is a good question. The last time we looked at it, the
C++ standard did not guarantee that a container of this type would
necessarily be associated with consecutive "packed" locations.
However, C and fortran arrays are always packed, so there is a
possibility that some implementations would be inconsistent.
In particular, it's not clear if you will always be able to
feed a C++ container to something like a legacy (or even "new" native C)
BLAS implementation. This may be what you were thinking about.

I don't think there is very much we can do about that.
No matter what you do, making an STL library implementation
talk to legacy codes will always require exposing some
detail of the STL implementation.

On the other hand, using STL containers for numerical
arrays is probably a bad idea anyway. That's why the
C++ standard includes valarray. Unfortunately, the
concensus seems to be that valarray is broken and
will never really work right (see discussions on the
mailing list at oonumerics.org). That leaves only
nonstandard implementations. But the state of the art
in nonstandard implementations is quite high. See
http://www.lsc.nd.edu/research/mtl/ for example.

In the end, you probably want to use STL for
one or more of the following reasons:

  1) standard, therefore portable and easy to use
  2) ready-made containers
  3) ready-made algorithms

If you use a third-party library, like MTL, you get
ready-made containers, so (2) is taken care of. If the
library containers are generic enough (MTL containers are),
then you get to use the STL algorithms that make sense.
That's the point of genericity, so you pretty much get (3)
for free too. That leaves (1), which is the sticky part.
With luck, genericity buys you a significant part
of (1) as well. But there will always the problem
of obtaining, building, installing, and tracking
for third-party products. Thankfully, it is not
too difficult when the authors show some
concern for these things.


Hope that answers your questions.
Thanks.

-- 
G. Jungman

Re: GSL

[Brian Gough]

Hi,

I am one of the developers.  Thanks for your questions, I will try to
answer them...

Jens Lenge writes:
 > ==> Is it the goal of GSL to build a "native C" equivalent of LAPACK
 > functionality?

We are trying to write a general scientific library which provides
similar functionality to existing fortran libraries.  However, we
don't intend to rewrite LAPACK because it is (a) huge and (b) of high
quality, so we could not improve on it.  We recommend using LAPACK for
large problems.  

We have provided some simpler routines for linear algebra, similar to
those found in Numerical Recipes.  These allow people to avoid
installing LAPACK if they only need to solve a small linear algebra
problem.

 > ==> Will it depend on (or include) an underlying BLAS library?

Yes.  We use blas internally in the library.

 > ==> If so, is that BLAS also native C?

Yes, we have one (partly written -- we still need to write a Level-3
BLAS).  We have a facility for substituting other BLAS libraries
(provided they conform to the Cblas standard).

 > ==> Who is actually doing GSL? Is Redhat involved?

There is a list of developers in the AUTHORS file of the distribution,
also accessible from the website http://sources.redhat.com/gsl/

Redhat is not directly involved in the development, but they provide
the hosting for the CVS access and website (much appreciated).

 > ==> When used as a basis for a C++ wrapper lib, will the GSL data structures
 > be compatible with STL containers or the C++ complex<double> class?

Our complex data type struct is { double d[2]; } where d[0] = real
part, d[1] = imag part.  This is compatible with most complex
datatypes in practice, e.g. {double re; double im}.  All the complex
number access is through macros so it is possible to recompile the
library with a different definition of complex if desired.

We have a design document gsl-design.texi which is in the doc/
directory of the gsl-0.7 tar file.

If you have any more questions feel free to ask.

regards
Brian Gough

Re: GSL

[Brian Gough]

Piero Foscari writes:
 > One thing isnt too clear, going for speed and Box-Muller is slowed down to 
 > stay reentrant... is then reentrance really a requirement along the whole 
 > project? I mean, other generators keep a state or even a whole set of 
 > values saved, why not do the same in this case?

There's a distinction between the functions in gsl_randist (variates
from distributions) and gsl_rng (generators).  The functions for
producing non-uniform variates don't maintain any state, since it
would tend to make them very complicated.

-- 
Brian Gough

Network Theory Ltd,
Commercial support for GSL --- http://www.network-theory.co.uk/gsl/

Re: GSL

[Andrea Riciputi]

Sorry for my previous post. Wrong address, wrong list.

Andrea.

On Thursday, Nov 13, 2003, at 12:39 Europe/Rome, Andrea Riciputi wrote:

> I've not installed 10.3 yet, but some time ago I submited a bug report 
> to Apple about a problem occurred during the testing phase of the 
> GSL's installation process. Apple aswered a couple of days ago saying 
> they can't reproduce the bug any more with the new compiler. I can 
> imagine that this mean GSL works well with Pather or at least I hope.
>
> Cheers,
> 	Andrea.
>

---
Andrea Riciputi

"Science is like sex: sometimes something useful comes out,
   but that is not the reason we are doing it" -- (Richard Feynman)

Re: GSL

[Andrea Riciputi]

I've not installed 10.3 yet, but some time ago I submited a bug report 
to Apple about a problem occurred during the testing phase of the GSL's 
installation process. Apple aswered a couple of days ago saying they 
can't reproduce the bug any more with the new compiler. I can imagine 
that this mean GSL works well with Pather or at least I hope.

Cheers,
	Andrea.

On Wednesday, Nov 12, 2003, at 05:00 Europe/Rome, Rodney Dyer wrote:

> Hey, has anyone tried the GNU Scientific Library on 10.3?  I'm getting 
> ready to upgrade and I'm interested in any problems that may be 
> involved.
>
> While I have your attention, any preferences for other libraries that 
> have LU and QR routines in C that do well with VERY large matrices?
>
> Thanks,
>
---
Andrea Riciputi

"Science is like sex: sometimes something useful comes out,
   but that is not the reason we are doing it" -- (Richard Feynman)

Re: GSL

[Brian Gough]

Andrew Royal writes:
 > DEar GSL programmers,
 > 
 > I would like to know if it is possible to use the GSL (originally
 > designed for Microsoft's Visual C++) for the free compiler devC++.  If
 > so, then how does one go about installing this.  I have looked
 > extensively on the web but can not find good documentation.  I am new to
 > programming and so I don't full understand all the jargon.  Are you able
 > to help?
 > 
 > In anticipation
 > 
 > Andrew

I don't use DevC++ myself so I have forwarded your message to the gsl
mailing list. You may want to ask on the Dev C++ mailing list if
nobody there knows.

regards
-- 
Brian Gough

Network Theory Ltd            
15 Royal Park                 
Bristol BS8 3AL               
United Kingdom                

Tel: +44 (0)117 3179309 
Fax: +44 (0)117 9048108              
Web: http://www.network-theory.co.uk/

Re: gsl

[Brian Gough]

Lotte Høj Andersen writes:
 > Hi Brian Gough,
 >  We are two students from Aalborg University, Denmark, and we have
 > a problem:
 >  we are doing a plastic finite element analysis, and we are using
 > the gsl-library. When we are doing our calculations the all vectors
 > and matrices that are allocated can not be free'd that is we are
 > using the command 'gsl_vector_free' but it makes no diffence.
 >  The calculations is performed in steps and we are using a number
 > of functions in each step, we are supposed to do about a 100 steps,
 > but since the allocated vectors and matrices are taking up a lot of
 > space (that is not beeing free'd) it is only possible to do about
 > 10 steps before the Harddisk is full.

Hi,

If the library is not working as described in the manual please can
you send a bug report to the mailing list
gsl-discuss@sources.redhat.com with an example program which
reproduces the problem. Thanks.

Brian Gough

Re: GSL

[Brian Gough]

John Aquilina writes:
 > i would like to install the GSL library in order to be able to
 > compute hypergeometric functions using the special fuinctions
 > package.
 >  What should i do to install?

The options for different platforms are listed on
http://sources.redhat.com/gsl

 > Do i have to download and extract all the library?

Yes, the library is installed as a single compiled package, not as
individual source files.

regards
Brian Gough

GSL

[John Aquilina]

i would like to install the GSL library in order to be able to compute
hypergeometric functions using the special fuinctions package.

What should i do to install? Do i have to download and extract all the
library?


thanks
John

 ---------------------------------
 John Aquilina
 Dept. of Mathematical Sciences
 University of Bath
 Claverton Down
 Bath BA2 7AY

 Office: 1W 3.10G
 Tel:    (+44) (0)1225 384738
 Mobile: (+44) (0)7761 380238
 E-Mail: J.Aquilina@bath.ac.uk
 ---------------------------------

Re: gsl

[Brian Gough]

Lotte Høj Andersen writes:
 > Hi Brian Gough,
 >  We are two students from Aalborg University, Denmark, and we have
 > a problem:
 >  we are doing a plastic finite element analysis, and we are using
 > the gsl-library. When we are doing our calculations the all vectors
 > and matrices that are allocated can not be free'd that is we are
 > using the command 'gsl_vector_free' but it makes no diffence.
 >  The calculations is performed in steps and we are using a number
 > of functions in each step, we are supposed to do about a 100 steps,
 > but since the allocated vectors and matrices are taking up a lot of
 > space (that is not beeing free'd) it is only possible to do about
 > 10 steps before the Harddisk is full.

Hi,

If the library is not working as described in the manual please can
you send a bug report to the mailing list
gsl-discuss@sources.redhat.com with an example program which
reproduces the problem. Thanks.

Brian Gough

Re: gsl

[Brian Gough]

Hiroyuki Kawakatsu writes:
 > i have some comments and questions regarding gsl (i am using
 > version 0.9.2 on windoze 2000).

Hi, 
Thanks for your comments.

 >  (1) the major problem with the current implementation is that
 > matrices are stored in row-major order. since i have been writing
 > codes that link with dll's built from fortran codes, i have assumed
 > all matrices will be stored in column-major order. for me to really
 > use gsl, i really need an option to store matrices in column-major
 > order.  for example, tnt (template numerical toolkit
 > http://math.nist.gov/tnt/ ) supports a fortran_matrix type in
 > addition to the row-major order c style matrix type and intel's
 > math kernel library has an argument (at least in cblas) to specify
 > the major order (isn't this the cblas standard's recommendation?)

Sorry, it's a C-library so we use arrays in the C convention for the
gsl_matrix object.

To support both conventions would either reduce the efficiency or
require separate handling of the two cases.  TNT uses C++ templates to
handle both cases separately, but this is not feasible in C.

The CBLAS standard does specify an argument for both conventions, and
the GSL CBLAS routines support that.

 > (2) in the documentation, it is mentioned that you can link with
 > atlas for cblas. is it possible to link with intel's math kernel
 > library's cblas? if so, it would be very nice if you could document
 > how to do so.

The section "Compiling and Linking" in Chapter 2 of the manual
describes how to link with an alternate CBLAS library.  Assuming the
library follows the CBLAS standard no changes are required to the
code, just to specify the library on the link options.

Intel's math library is not free software, so we don't comment on it
specifically.  We recommend ATLAS, which is free software.

 >  (3) i am a bit disappointed with the optimization routines which
 > seem to be fairly incomplete at this stage. (i am quite happy with
 > the pseudo-random number generator routines which are up-to-date.)
 > i tried to test the non-linear least squares routine using the nist
 > problems http://www.itl.nist.gov/div898/strd/nls/nls_main.shtml but
 > don't seem to get any accuracy at all. i am skeptical with the
 > claim in the documentation that the routines are based on minpack
 > since the dll version of minpack passes most of the nist test
 > problems. of course, there may be some problems with my c++ code
 > (which i had to convert to row-major ordering). if you have tested
 > the code with the nist problems, i would like to see the results
 > (and code).

I haven't run the NIST problems, but you'll find a working test
program in the multifit/ directory.

If you still have problems after checking the program send the code
the gsl-discuss@sources.redhat.com list for comments.

 >  (4) for the percentile function
 > gsl_stats_quantile_from_sorted_data(), i suggest that you provide
 > options as described in
 >  Rob J. Hyndman and Yanan Fan (1996) "Sample quantiles in
 > statistical packages," American Statistician, 50 361-365.

I'll take a look at that paper next time I'm at the library.

best regards,

Brian Gough

Re: gsl

[Edwin Robert Tisdale]

Brian Gough wrote:

> Hiroyuki Kawakatsu writes:

> >  (1) the major problem with the current implementation is that
> > matrices are stored in row-major order.
> > Since I have been writing codes that link with dll's built from
fortran codes,
> > I have assumed that all matrices will be stored in column-major
order.
> > For me to really use gsl, I really need an option
> > to store matrices in column-major order.
> > For example, the Template Numerical Toolkit (TNT)
> >
> >     http://math.nist.gov/tnt/
> >
> > supports a fortran_matrix type
> > in addition to the row-major order c style matrix type
> > and Intel's math kernel library has an argument (at least in cblas)
> > to specify the major order
> > (isn't this the cblas standard's recommendation?)
> >
> > Sorry, it's a C-library
> > so we use arrays in the C convention for the gsl_matrix object.
> >
> > To support both conventions would either reduce the efficiency
> > or require separate handling of the two cases.
> > TNT uses C++ templates to handle both cases separately
> > but this is not feasible in C.
> >
> > The CBLAS standard does specify an argument for both conventions
> > and the GSL CBLAS routines support that.
> >
> > (2) in the documentation, it is mentioned that
> > you can link with atlas for cblas.
> > Is it possible to link with intel's math kernel library's cblas?
> > If so, it would be very nice if you could document how to do so.
> >
> > The section "Compiling and Linking" in Chapter 2 of the manual
> > describes how to link with an alternate CBLAS library.
> > Assuming that the library follows the CBLAS standard,
> > no changes are required to the code
> > just to specify the library on the link options.
> >
> > Intel's math library is not free software,
> > so we don't comment on it specifically.
> > We recommend ATLAS, which is free software.

Hi Hiroyuki,

Neither your computer nor the C computer programming language
know anything about rows or columns.
Rows and columns exist only in the mind of the computer programmer.
The only reason that you can pass two dimensional arrays
from C programs to Fortran subprograms or
from Fortran programs to C functions is that
most Fortran and C compilers store two dimensional arrays
in exactly the same way -- namely as arrays of one dimensional arrays.
The only difference between multidimensional arrays in C and
multidimensional arrays in Fortran (besides the default base index)
is that Fortran subscripts appear in order from left to right and
C subscripts appear in reverse order from right to left.
Unfortunately, both Fortran and C programmers associate
the leftmost subscript of a two dimensional array
with the row index of a matrix and associate
the rightmost subscript of a two dimensional array
with the column index of a matrix so that
the same matrix object implemented as a two dimensional array
with one language appears to be its transpose in the other language.

What this means is that you never need to transpose a matrix
to call a Fortran subprogram from C (or to call a C function from
Fortran.)
All you need to do is to call the correct function with the correct
arguments.
For example, if you wish to compute

    alpha*A*B + beta*C --> C

where C in an m by n matrix,
you would first transpose the expression

    alpha*B^T*A^T + beta*C^T --> C^T

then transpose each matrix in the expression

    alpha*B*A + beta*C --> C.

In other words, you would still call xGEMM
but swap A and B ( and m and n) in the argument list.

Re: gsl

[Brian Gough]

Hiroyuki Kawakatsu writes:
 > i tried to test the non-linear least squares routine using the nist
 > problems http://www.itl.nist.gov/div898/strd/nls/nls_main.shtml but
 > don't seem to get any accuracy at all.

You're right, I've run a couple of the nist tests and there is a bug
which prevents convergence for some problems.  Thanks. I will fix that.