Re: gsl

[Edwin Robert Tisdale <E dot Robert dot Tisdale at jpl dot nasa dot gov>]

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.