Re: dgemm subroutine in BLAS - I think I've cracked the difference, please confirm

From: glen herrmannsfeldt (gah_at_ugcs.caltech.edu)
Date: 02/14/05 

Date: Sun, 13 Feb 2005 15:18:38 -0800

Dawn Minnis wrote:

> Oh do you know I may have just got a blow to the head there. By any chance
> does LDA,B,C define the size of one dimension of the BIG BIG BIG outer
> matrix and m,n,k define the dimensions of the actual matrix to be worked
> upon. And going on that, LDA, B, C explain what the jump is between one
> element in row eg 3 and the next element (ie in the next col)?

> Ok, if thats right then I get it. I was actually working on it the wrong
> way round thinking that m,n,k defined the outer matrix and LDA,B,C defined
> the inner size.

> So I build a matrix first of all, being of size A=(LDA, k) B=(LDB, n) and
> C=(LDC, n) and then only work with the submatrix A=(m, k) B=(k, n) and C=(m,
> n).

> ok, (as you can see I'm talking this out as I understand it so even I can
> see what I'm talking about), so in the current state of memory allocation,
> given that there is a function malloc() which I am using to allocate the
> required memory for A,B and C, is there acually any need to define LDA,B,C
> anymore? I mean, if you are wanting to perform multiplication of matrix A*B
> then why bother defining a larger matrix to begin with that takes up
> memory??

Well, that is a completely different question.

Now that you know all about Fortran arrays and memory use, everything
is completely different in C.

First, C stores arrays with the rightmost subscript varying fastest,
instead of leftmost as in Fortran. One reason is that it is closer
to mathematical conventions, and the other has to do with the connection
between arrays and pointers in C.

While newer Fortran versions actually supply dymamic allocation
for arrays, C supplies the tools but leaves out much of what you
need for this problem.

C does supply statically allocated multidimensional arrays, though
some would disagree with calling them that. For dynamic allocation
malloc() can only allocate one in one dimension.

Using C's static array allocation, you can say something like

double a[100][100],b[100][100],c[100][100];

which will declare three two-dimensional arrays somewhat
like in Fortran, except with the dimensions in the opposite
order, and only constants (or constant expressions) are allowed
for the dimensions. In this case, all you learned last week will
be extra useful. With a little luck you can call Fortran subroutines
using such arrays.

If you want to use malloc, things get more complicated.
You can, for example, do:

double *a, *b, *c,alpha=2,beta=0;
char no='n';
scanf("%d%d%d",&k,&l,&m);
/* mk, kn, mn */

a=malloc(m*k*sizeof(*a));
b=malloc(k*n*sizeof(*b));
c=malloc(n*m*sizeof(*c));

/** now you can store in these arrays in the same order as
Fortran would: **/

for(i=0;i<k;i++) for(j=0;j<m;j++) {
    scanf("%lf",a+i*m+j);
    }
for(i=0;i<n;i++) for(j=0;j<k;j++) {
    scanf("%lf",a+i*k+j);
    }

/** the exact way to call a Fortran subroutine from C is implementation
dependent. This is one possibility. Character strings are especially
questionable, and unfortunately they are needed here. **/

DGEMM ( &no, &no, &m, &n, &k, &alpha, a, &m, B, &k, &beta, &c,&m);

for(i=0;i<m;i++) {
    for(j=0;j<n;j++) printf("%20.10g ",c[j*m+i]);
    printf("\n");
    }

Note that C is no help in computing the position in the
array, in the way that Fortran does it.

-- glen

Relevant Pages

  • Re: INTERFACE problem
    ... I have an integer at Fortran side, pass it by reference to C function, treat here as pointer, allocate memory and this C-pointer at Fortran side now is an integer which has represents address allocated at C side. ... I did something like this many years ago to allow dynamic allocation ... Then it was converted into a subroutine with the arrays and ... it called the sub by way of Pascal with those "pointers". ...
    (comp.lang.fortran)
  • dgemm subroutine in BLAS - I think Ive cracked the difference, please confirm
    ... see what I'm talking about), so in the current state of memory allocation, ... but I have a problem when it comes to applying the dimensions. ... > In the Fortran 66 days it was because dynamic allocation didn't exist, ... > you allocate arrays big enough for the largest problem you might run into. ...
    (comp.lang.fortran)
  • Re: Speed penalty for using allocatable arrays?
    ... "Leading dimension" means dimension 1, ... Fortran compilers, rather than having being left those writing Fortran ... allocation available in Fortran 90/95. ... the arrays are allocated once at the ...
    (comp.lang.fortran)
  • Re: Fortran memory allocation (stack/heap) issues
    ... > rather than Fortran, ... dynamic allocation, and relatively little stack allocation. ... value return and arrays by reference. ...
    (comp.lang.fortran)
  • Fortran and Python (was Re: Statement function host association)
    ... Are trying to sneak Python into Fortran? ... procedures taking allocatable arrays as arguments and returned ... automatic allocation of arrays in an expression ...
    (comp.lang.fortran)