Multiplying a subset of matrix to another matrix using lapack/blas

Question

I want to multiply a subset of matrix A to another matrix, using dgemm or anyother lapack/blas function. I think that since elements of the submatrix may not be contiguous, I cannot use dgemm directly without copying the submatrix to another space. So, when this submatrix is itself large, it may be very inefficient to the extent that I think it may be better for me to write the code for multiplication for this specific problem in C. Since copying and then using the lapack/blas itself, may not be efficient at all. I am using lapack/blas in matlab as a mex file.

My questions are

1- Is there any function of lapack/blas which can work on submatrices in multiplication? 2- If it isn't, is it better to write the code for multiplication directly or is it better to copy submatrix to another matrix and use dgemm?

ztik · Accepted Answer · 2017-09-26T08:37:05.247

Actually dgemm is designed for submatrix multiplication. You just need to make right use of the starting pointers for each matrix and the arguments LDA, LDB, LDC.

The C variant of BLAS is :

void cblas_dgemm (const CBLAS_LAYOUT layout, const CBLAS_TRANSPOSE TransA, const CBLAS_TRANSPOSE TransB, const int M, const int N, const int K, const double alpha, const double *A, const int lda, const double *B, const int ldb, const double beta, double *C, const int ldc);

Say that you have matrices:

A(15x10)
B(10x20)
C(15x20)

Calling dgemm for Column Major matrix storage :

cblas_dgemm (CblasColMajor, CblasNoTrans, CblasNoTrans, 15, 20, 10, 1., A, 15, B, 10, 1., C, 15);

Suppose that you need to call dgemm passing the submatrices:

As(3x2) starting at point (2,1) of A
Bs(2x5) starting at point (3,5) of B
Cs(3x5) starting at point (4,2) of C

The N, M, K will change to 3, 5, 2, but the LDXs will remain same as above. Then you have to pass the right pointers to dgemm so that they will point to the start of each submatrix. Since you have C numbering you have to subtract one from each coordinate.

As starting point is A + (1+0*15)
Bs starting point is B + (2+4*10)

Cs starting point is C + (3+1*15)

cblas_dgemm (CblasColMajor, CblasNoTrans, CblasNoTrans, 3, 5, 2, 1., A+1, 15, B+42, 10, 1., C+18, 15);

The idea of N LDA is to say that I have a matrix A(LDA,*) but I will use the upper submatrix As(N,*). In the examples case you do not want to use the upper submatrix but some other inside A. In this case you create a new pointer A+1 to matrix. Now As is an upper submatrix of A+1.

Similarly calling the Fortran's original dgemm function from C will be

    char NoTrans = `N`;
    int N = 3;
    int M = 5;
    int K = 2;
    int LDA = 15;
    int LDB = 10;
    int LDC = 15;
    double alpha = 1.0;
    double beta = 1.0;
    dgemm (&NoTrans, &NoTrans, N, M, K, alpha, A+1, LDA, B+42, LDB, beta, C+18, LDC);

Thank you very much. Could you please explain a little more what it means to use it this way? How blas_dgemm finds the next elements of the submatrix given that it doesnot know the start of the matrix and size of each column? and can I use this in lapack/blas in matlab? sorry I am so new to this. — user85361, Sep 26 '17 at 08:13
@user85361 I just mentioned cblas equivalent because you mentioned `C`. You can use original Fortran's `dgemm` in your mex file as well. I will update the answer — ztik, Sep 26 '17 at 08:29
How to multiply two submatrices when we have [0/1] indices for row and columns, so they may not be contiguous? — user85361, Oct 31 '17 at 02:10
In matlab you can have logical indexing, e.g., lind = myVec>0. In this way you have indices which 1 means it is in the computation and 0 means it's not. — user85361, Oct 31 '17 at 21:59

Multiplying a subset of matrix to another matrix using lapack/blas

1 Answers1