Re: simplicity vs. efficiency of algorithms in GSL

[Gerard Jungman <jungman@lanl.gov>]

On Fri, 2008-10-10 at 22:13 +0100, Brian Gough wrote:
> > 
> > It's a draft standard for a C LAPACK interface.
> 
> Interesting, but it's column major only.

Yup.

Here's the way I see it. Anybody who uses lapack already
faces this problem; if you want to use lapack, you have
to get your ducks lined up in columns at the start.
People must already be doing this. So it makes sense to
me that a library (GSL) should encapsulate whatever they
are doing, so they don't have to do it over and over again.

This is the problem that we face. The C/C++ community must
decide for themselves what it means to "use lapack". Does it
mean that column-major is required? If so, how does
one manage both the row-major and column-major data
types? A library should presumably provide both.
What are the relative costs of runtime conversion?
For many uses that could be a viable option, if
people are forced into using both types in code.
People who use only lapack and blas functions
probably don't care and don't need to know how
their data is laid out, and could be given the
column-major type by default.

How should a library support both types? Clearly,
the data layer need not change, only the access
patterns need to be abstracted. Those patterns
can be parametrized in such a way as to make
code generation easy. Personally, I would
prefer C++ templates for this, but we can
use whatever C code-generation kludge people
happen to fancy.

Are these problems best solved as part of a general
"interface to fortran", or are they specific to
matrix operations? Is any of this addressed in
the new fortran standard (fortran-2003 apparently
has specific standardization for inter-language
support, including stuff like passing structs
back and forth, but I don't know the details).
Obviously they can't tell you how to lay out
your data, but any guidance in this area could
be helpful.

Finally, we should contact the people working on
this interface standard. It may be just that one
guy. He might appreciate some discussion.

As people are probably aware, the cblas interface
is data-order agnostic; you can treat any matrix
as row-major or column-major, using a flag. It is
unfortunate that the proposed lapack standard does
not adopt this idea. But we may have to accept the
fact that lapack is essentially legacy code,
and that's just the way it is.

Let's talk to them and find out.

--
G. Jungman