Search code examples
matlabmatrixmatrix-multiplicationlapackblas

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

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?

Solution

  • 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);