How Does cblas_dgemm Handle Matrix Operations in C?

[KStolen]

Hi guys, I'm having trouble understanding how this routine works.

cblas_dgemm is a BLAS function that gives C <= alpha*AB + beta*C
where A,B,C are matrices and alpha, beta are scalars.

 void cblas_xgemm (
     const enum CBLAS_ORDER Order,
     const enum CBLAS_TRANSPOSE TransA,
     const enum CBLAS_TRANSPOSE TransB,
     const int M,
     const int N,
     const int K,
     const SCALAR alpha,
     const TYPE * A,
     const int lda,
     const TYPE * B,
     const int ldb,
     const SCALAR beta,
     TYPE * C,
     const int ldc)

It takes 14 parameters, listed http://www.psatellite.com/matrixlib/api/lapack.html" .

http://en.wikipedia.org/wiki/General_Matrix_Multiply" on wikipedia.

I don't understand what the "major stride" (lda,ldb,ldc) is or how it works, despite the explanations given on both sites.

Here are my example matrices:

A = \begin{bmatrix}1 <br /> &amp;1 &amp;1 &amp;1 \\ 1 <br /> &amp;1 &amp;1 &amp;1 \\ 1<br /> &amp;1 &amp;1 &amp;1 \\ 1<br /> &amp;1 &amp;1 &amp;1 <br /> \end{bmatrix}<br />

<br /> B = \begin{bmatrix}1 <br /> &amp;1 &amp;1 &amp;1 \\ 1 <br /> &amp;1 &amp;1 &amp;1 \\ 1<br /> &amp;1 &amp;1 &amp;1 \\ 1<br /> &amp;1 &amp;1 &amp;1 <br /> \end{bmatrix}<br />

<br /> C = \begin{bmatrix}0 <br /> &amp;0 &amp;0 &amp;0 \\ 0 <br /> &amp;0 &amp;0 &amp;0 \\ 0<br /> &amp;0 &amp;0 &amp;0 \\ 0<br /> &amp;0 &amp;0 &amp;0 <br /> \end{bmatrix}<br />

I want to be able to take the bottom quarter submatrices of A and B, multiply them and add them to C such that

<br /> C = \begin{bmatrix}0 <br /> &amp;0 &amp;0 &amp;0 \\ 0 <br /> &amp;0 &amp;0 &amp;0 \\ 0<br /> &amp;0 &amp;2 &amp;2 \\ 0<br /> &amp;0 &amp;2 &amp;2 <br /> \end{bmatrix}<br />

Here's an excerpt from my code:

//where n is the matrix size, in this case 4
void Multiply(int n, int blockSize, double** a, double** b, double** c)
{
       cblas_dgemm(CblasRowMajor,CblasNoTrans, CblasNoTrans, n, n , n , 1.0, a[0], n, b[0], n, 1.0, c[0], n);
}

This code successfully multiplies A by B and gets C, a matrix filled with 4's.

I'm thinking something like:

//where n is the matrix size, in this case 4
void Multiply(int n, int blockSize, double** a, double** b, double** c)
{
       cblas_dgemm(CblasRowMajor,CblasNoTrans, CblasNoTrans, 2, 2 , 2 , 1.0, a[1], n, b[1], n, 1.0, c[1], n);
}

but that returns

<br /> C = \begin{bmatrix}0 <br /> &amp;0 &amp;0 &amp;0 \\ 2 <br /> &amp;2 &amp;2 &amp;2 \\ 0<br /> &amp;0 &amp;0 &amp;0 \\ 0<br /> &amp;0 &amp;0 &amp;0 <br /> \end{bmatrix}<br />

instead.

I think it's because I don't understand the lda, ldb and ldc parameters. Is it possible to get what I want here, through only using cblas_dgemm? Obviously, I could iterate over the arrays with loops but I'd prefer not to have to do that.

 

[KStolen]

Well I realized what I was doing wrong, so I thought I'd put the description here for anyone who needs to know.

lda, ldb and ldc (the strides) are not relevant to my problem after all, but here's an explanation of them :

The elements of a matrix (i.e a 2D array) are stored contiguously in memory. However, they may be stored in either column-major or row-major fashion. The stride represents the distance in memory between elements in adjacent rows (if row-major) or in adjacent columns (if column-major). This means that the stride is usually equal to the number of rows/columns in the matrix.

Matrix A =
[1 2 3]
[4 5 6]
Row-major stores values as {1,2,3,4,5,6}
Stride here is 3

Col-major stores values as {1, 4, 2, 5, 3, 6}
Stride here is 2


Matrix B =
[1 2 3]
[4 5 6]
[7 8 9]

Col-major storage is {1, 4, 7, 2, 5, 8, 3, 6, 9}
Stride here is 3


However, my issue was very simple in the end. I just didn't understand pointers. When you reference an array in C, what you are really referencing is the memory location of the first element of that array. Hence, to perform the operation I required, simply do this:

//where n is the matrix size, in this case 4
//if blocking, replace M,N,K with blockSize
void Multiply(int n, int blockSize, double** a, double** b, double** c)
{
       cblas_dgemm(CblasRowMajor,CblasNoTrans, CblasNoTrans, 2, 2 , 2 , 1.0, a[2]+2, n, b[2]+2, n, 1.0, c[2]+2, n);
}

 

[Mark44]

The elements of a matrix (i.e a 2D array) are stored contiguously in memory. However, they may be stored in either column-major or row-major fashion. The stride represents the distance in memory between elements in adjacent rows (if row-major) or in adjacent columns (if column-major). This means that the stride is usually equal to the number of rows/columns in the matrix.

I'm reasonably sure that the elements of a two-dimension array are stored in row-major order in C and other C-based languages. If I remember correctly, Fortran does things differently, storing elements of a matrix in column-major order.

 

[KStolen]

That's right Mark. Because BLAS is written in Fortran, matrix row/column order is something you should be aware of when passing it matrices from C.