Fri Mar 31 23:19:00 GMT 2006
Gerard Jungman writes:
> - We discussed creating a lapack interface, similar to the
> blas interface. The basic limitation was lack of manpower
> and some statements about the difficulty of dealing with
> workspace allocation/management, which I did not understood
> at the time, and still don't understand.
The basic problem with a LAPACK interface is there is no clean way to
do it. I've tried this several times. As we remember, with BLAS you
have the 'CBLAS trick' to get the row-major versions of Fortran BLAS
by using the BLAS transpose argument, which is conveniently provided
in all necessary functions.
The LAPACK functions typically don't have a transpose argument and are
hard-coded for column-major. This means either you have to use
column-major in C (which is inconvenient if you have other C-style
matrices around) or you need to provide transposed matrices, either by
logically solving the transposed system or physically transposing them
before each call and then transposing them back (which is what Numpy
If the matrix has M=N or N=TDA then you can do the transpose in-place,
so you don't need extra memory. But in the general case M!=N, N!=TDA
you have to allocate a copy of the matrix. Obviously for large
matrices this is wasteful of memory, but these are the sorts of
systems you would want to solve.
Already you are looking at 3 possible types of wrapper function for
each LAPACK function to handle the different cases. But it's more
complicated than that because for functions with multiple arguments
you might have different cases for each argument.
So either there has to be some radical simplifying assumption like
"always allocate memory and transpose" or "only handle matrices with
N=TDA" or "make the user pass everything as column-major". But people
will surely complain whichever is chosen.
Until LAPACK has transpose arguments to make a 'LAPACK trick' possible
it's a lost cause and is probably easier for people call LAPACK
directly on a case-by-case basis.
More information about the Gsl-discuss