Re: LAPACK vs. GSL matrix algebra

[Linas Vepstas <linas at austin dot ibm dot com>]

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