LAPACK vs. GSL matrix algebra

LAPACK vs. GSL matrix algebra

[Linas Vepstas]

Hi,

What is the current thinking with adding to the matrix functionality 
in GSL?   I need to find the eigenvectors of a non-symmetric real
matrix, which is something that GSL doesn't currently support.

My choices seem to be
-- try to remember how to code in fortran
-- try to implement something in C, with an eye towards submitting 
   it to GSL.

Any suggestions, words of encouragement to go one way or the other?


--linas

Re: LAPACK vs. GSL matrix algebra

[kaska]

you can easily access the fortran routines from C.
i always write C routines that access the fortran
functions from LAPACK that i need.  as long as you
keep in mind that fortran stores matrices column-wise
and C does it row-wise.  but most fortran routines
allow you to set the 'transpose' option so you can
easily pass in your C matrix.  
that way you don't have to write your own matrix
routines.  the fortran lapack is pretty good and it
has been extensively tested over the years.
if you need example files, let me know.
k



On Mon, 27 Sep 2004, Linas Vepstas wrote:

> 
> 
> Hi,
> 
> What is the current thinking with adding to the matrix functionality 
> in GSL?   I need to find the eigenvectors of a non-symmetric real
> matrix, which is something that GSL doesn't currently support.
> 
> My choices seem to be
> -- try to remember how to code in fortran
> -- try to implement something in C, with an eye towards submitting 
>    it to GSL.
> 
> Any suggestions, words of encouragement to go one way or the other?
> 
> 
> --linas
> 

Re: LAPACK vs. GSL matrix algebra

[C J Kenneth Tan -- OptimaNumerics]

Linas,

You don't have to remember how to code in Fortran to use LAPACK. 

BLAS and LAPACK are available for most platforms these days.  If you
use BLAS and LAPACK, you can always use the free versions from Netlib.
When you need speed, you can switch to optimized versions of BLAS and
LAPACK.  High performance implementations of these libraries can
deliver very significant speedup.  We have seen improvements of
hundreds to thousands of percents depending on the routine and
architecture, compared to hardware vendors' code!


Kenneth Tan
-----------------------------------------------------------------------
C. J. Kenneth Tan, Ph.D.
OptimaNumerics Ltd.
E-mail: cjtan@OptimaNumerics.com      Telephone: +44 798 941 7838    
Web: http://www.OptimaNumerics.com    Facsimile: +44 289 066 3015 
-----------------------------------------------------------------------
This e-mail (and any attachments) is confidential and privileged.  It
is intended only for the addressee(s) stated above.  If you are not an
addressee, please accept my apologies and please do not use,
disseminate, disclose, copy, publish or distribute information in this
e-mail nor take any action through knowledge of its contents: to do so
is strictly prohibited and may be unlawful.  Please inform me that
this e-mail has gone astray, and delete this e-mail from your system.
Thank you for your co-operation.
-----------------------------------------------------------------------




On 2004-09-27 11:20 -0500 Linas Vepstas (linas@austin.ibm.com) wrote:

> Date: Mon, 27 Sep 2004 11:20:18 -0500
> From: Linas Vepstas <linas@austin.ibm.com>
> To: gsl-discuss@sources.redhat.com
> Subject: LAPACK vs. GSL matrix algebra
> 
> 
> 
> Hi,
> 
> What is the current thinking with adding to the matrix functionality 
> in GSL?   I need to find the eigenvectors of a non-symmetric real
> matrix, which is something that GSL doesn't currently support.
> 
> My choices seem to be
> -- try to remember how to code in fortran
> -- try to implement something in C, with an eye towards submitting 
>    it to GSL.
> 
> Any suggestions, words of encouragement to go one way or the other?
> 
> 
> --linas
> 
> 
> 

Re: LAPACK vs. GSL matrix algebra

[Gerard Jungman]

On Mon, 2004-09-27 at 10:35, kaska@yatsen.cas.McMaster.CA wrote:
> you can easily access the fortran routines from C.
> i always write C routines that access the fortran
> functions from LAPACK that i need.

GSL should provide an interface to these routines,
the way it does for BLAS. We thought about it years
ago, but it never got done. There were some issues,
but probably no show-stoppers.

This may be the most important open issue in GSL.
The question comes up quite often here, and the
usual answer, "well, you'll have to figure out how
to connect to LAPACK yourself", gets embarrassing.

By the way, I would appreciate any comments from
people using an external BLAS with the GSL interfaces.
In principle it should be as simple as replacing
"-lgslcblas" with "-lmycblas". But I have never
heard any complaints or discussion of any kind,
which makes me think that nobody is doing
it this way.


-- 
Gerard Jungman <jungman@lanl.gov>
Los Alamos National Laboratory

Re: LAPACK vs. GSL matrix algebra

[Linas Vepstas]

On Mon, Sep 27, 2004 at 12:06:03PM -0600, Gerard Jungman was heard to remark:
> On Mon, 2004-09-27 at 10:35, kaska@yatsen.cas.McMaster.CA wrote:
> > you can easily access the fortran routines from C.

Indeed.  The hardest part was reading the lapack doco's
and creating  the C prototype for the needed routine.

The other hard part was figuring out the libraries to link to:

.o:
   cc -o $@ $< -llapack -lf77blas -latlas -lg2c -lm

.c.o:
   cc -g -O2 -c $<


> GSL should provide an interface to these routines,
> the way it does for BLAS. We thought about it years

Yes, it should.  The C prorotype for fortran routines
is kind-of ugly.

Did you auto-gen the cblas code and/or headers, or 
are these hand-written?

--linas

Re: LAPACK vs. GSL matrix algebra

[Gerard Jungman]

On Mon, 2004-09-27 at 14:02, Linas Vepstas wrote:
> 
> Did you auto-gen the cblas code and/or headers, or 
> are these hand-written?

In short, it was all done by hand.

The implementation of the GSL native cblas/ was extremely
tedious, but obviously we don't need that step for getting
lapack interfaces.

The stuff in the blas/ directory was not hard, and since
cblas was a stationary target it seemed reasonable to
just do it all by hand. It think it took me one long
afternoon to get it all in, plus some later cleanup.

Questions:

 Is there a C-LAPACK standard, in the same way that
 there is a C-BLAS standard? I think there was some
 work done, but I don't know if it was ever finished.

 Would GSL target a C-LAPACK standard or the standard
 fortran? It is probably more work to target fortran.

 What is the usage model? Suppose somebody wants to
 use GSL interfaces over an external lapack. The external
 lapack is presumably implemented in terms of an external
 blas, probably both in fortran. What if the user also
 uses explicit blas functions (also through a GSL interface)?
 How do we guarantee that the blas called by lapack functions
 is consistent with the blas available through GSL interfaces.
 The goal is to write code using GSL interfaces and be able
 to switch in any conformant linear algebra sub-system with
 at most a re-compile, and perhaps only with a re-link.

 Hasn't somebody done all this already??? And can't
 we arrange so that all the work is done through ATLAS?


We've had these discussions here before. I wish there
was some actual plan to make this all work right.

-- 
Gerard Jungman <jungman@lanl.gov>
Los Alamos National Laboratory

Re: LAPACK vs. GSL matrix algebra

[Brian Gough]

Gerard Jungman writes:
 > By the way, I would appreciate any comments from
 > people using an external BLAS with the GSL interfaces.
 > In principle it should be as simple as replacing
 > "-lgslcblas" with "-lmycblas". But I have never
 > heard any complaints or discussion of any kind,
 > which makes me think that nobody is doing
 > it this way.

In the past I run have the test suite with ATLAS as well as GSL CBLAS
and it worked fine.

-- 
Brian Gough

Re: LAPACK vs. GSL matrix algebra

[Linas Vepstas]

On Mon, Sep 27, 2004 at 02:54:50PM -0600, Gerard Jungman was heard to remark:
> On Mon, 2004-09-27 at 14:02, Linas Vepstas wrote:
> > 
> > Did you auto-gen the cblas code and/or headers, or 
> > are these hand-written?
> 
> In short, it was all done by hand.

Hmm. Trying to figure out if I should volunteer to do this.

>  Is there a C-LAPACK standard, in the same way that
>  there is a C-BLAS standard? I think there was some
>  work done, but I don't know if it was ever finished.

Google only seems to hit the f2c'ed-up version of lapack.
Whatever work was done towards that standard seems to be burried
somewhere.

The f2c'ed c-lapack binding seem to be (superfically?) identical 
to the fortran bindings; i.e. fairly ugly looking to the 
c-programmers eye.

I could try to hand-build a more acceptable binding, but would
like to have active review for that.

>  Would GSL target a C-LAPACK standard or the standard
>  fortran? It is probably more work to target fortran.

More work in what way?  The bindings from above look to be
the same, more or less; any 'hard work' would be to 
swap out blas from underneath, if that was desired.

>  What is the usage model? Suppose somebody wants to
>  use GSL interfaces over an external lapack. The external
>  lapack is presumably implemented in terms of an external
>  blas, probably both in fortran. 

Bing.  So I guess that you are thinking of just literally 
copying c-lapack into the gsl cvs tree, and using that as 
the starting point?

Are there any license problems? lapack seems to be either 
bsd licnese, or public domain; I can't tell which; gsl is gpl'ed.

>  What if the user also
>  uses explicit blas functions (also through a GSL interface)?
>  How do we guarantee that the blas called by lapack functions
>  is consistent with the blas available through GSL interfaces.
>  The goal is to write code using GSL interfaces and be able
>  to switch in any conformant linear algebra sub-system with
>  at most a re-compile, and perhaps only with a re-link.

The "professional" way to solve this is to do it at run-time, 
by calling through a table of function pointers.  By default,
the table points to gsl-cblas, but either the whole teable or
individual routines in it can be replaced by the user, as desired.

>  Hasn't somebody done all this already??? And can't
>  we arrange so that all the work is done through ATLAS?

Or maybe the function pointer table would be initialized 
with atlas routines I guess ... 

Do you want a sketch of what the function-pointer table might 
look like?

--linas

Re: LAPACK vs. GSL matrix algebra

[Gerard Jungman]

On Wed, 2004-09-29 at 11:54, Linas Vepstas wrote:
> 
> Hmm. Trying to figure out if I should volunteer to do this.

You would be our hero.

> Google only seems to hit the f2c'ed-up version of lapack.
> Whatever work was done towards that standard seems to be burried
> somewhere.

The interfaces expressed in the ATLAS clapack.h might be
(part of) what I am thinking of. Clint Whaley must have some
notion of what to do for a complete wrapper, although
the ATLAS clapack.h interfaces have only some of the
functions. I vaguely remember some discussion about
standardization of C interfaces years ago. The C-BLAS
interface draft came out around 1998, and there may
have been some discussion of lapack interfaces at that
time.

> I could try to hand-build a more acceptable binding, but would
> like to have active review for that.

Right. The logical thing would be to create a generic C
interface and then target GSL to that. That means you
would effectively become the C-LAPACK standardizer.
We would need to get some input on that, probably at
least from Whaley if nobody else.

I guess the main issue is handling the workspace stuff,
although I haven't thought about it. It needs to be simple
but flexible. Brian thought about this about 5 years ago,
so he should say something...

> Bing.  So I guess that you are thinking of just literally 
> copying c-lapack into the gsl cvs tree, and using that as 
> the starting point?

It is probably necessary to have a "native" GSL implementation,
in the same way that we have a native cblas, just so people
have something that works out of the box. I originally wanted
to just dump a fortran blas implementation into a subdirectory
of GSL and use the cblas reference implementation to provide
a cblas; nasty and inefficient, but hardly any worse than our
hand-rolled cblas. Brian argued that we should roll our own,
I think based on the grodiness of having to carry around
and build fortran sources as an integral part of GSL.

Copying clapack into the source tree is one way to get the
functions to build out of the box. But I'm open to any ideas.
Also, it doesn't solve the general problem of creating
reasonable C interfaces; we would probably have to do
that anyway.

I wish that ATLAS had full lapack functionality. Then at
least we would have that to look at. I don't know what
his plans are in this regard.

> Are there any license problems? lapack seems to be either 
> bsd licnese, or public domain; I can't tell which; gsl is gpl'ed.

lapack itself is certainly free enough. The rpm that I have was
built by Redhat, from the netlib sources, and the "License" entry
says "Freely distributable". I couldn't find any markings on
the f2c'ed clapack sources, but I'm guessing that they can't
be more restrictive than lapack and f2c themselves.

> The "professional" way to solve this is to do it at run-time, 
> by calling through a table of function pointers.  By default,
> the table points to gsl-cblas, but either the whole teable or
> individual routines in it can be replaced by the user, as desired.

I guess that makes sense. It makes me a little nervous.

> Or maybe the function pointer table would be initialized 
> with atlas routines I guess ... 

In principle the libraries all use the standard cblas_XXX names.
How do you get all these linked together so that they are
available together at run time? Is this a dumb question?
This is the main reason that I thought the link-time
solution was natural.

> Do you want a sketch of what the function-pointer table might 
> look like?

Sure. It needs to be thread-safe, but that's probably not hard.
How do you solve the shared name problem?

-- 
Gerard Jungman <jungman@lanl.gov>
Los Alamos National Laboratory

Re: LAPACK vs. GSL matrix algebra

[Brian Gough]

As a test case, see if you can wrap DGESVX (AX=B) in a meaningful way.
I think it always ends up as "XA=B" in C, regardless of TRANS='N' or
'T'.

It looks like Clint Whaley made his own implementation of the small
number of LAPACK routines in ATLAS so that they could have an
additional CBLAS_ORDER argument.

-- 
Brian Gough

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