Can MEX BLAS library be used for native double matrix in C?

I'm writing C MEX file to do vector-matrix multiplication,

clear
clc

mex simpledgemv.c -R2017b -lmwblas
mex simpledgemv_native.c -R2017b -lmwblas

A = [1 0;0 1;2 0];
b = [3;4];
sol = A*b;
sol_blas = simpledgemv(A, b);    % Works well, 3.000000 4.000000 6.000000
sol_blasnative = simpledgemv_native();  % 0.000000 0.000000 0.000000

where, simpledgemv accepts arguments from MATLAB

#include "mex.h"
#include "matrix.h"
#include "blas.h"

#if !defined(_WIN32)
#define dgemv dgemv_
#endif

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    char *BLAS_NOTRANS = "N";
    const double dalpha = 1.0, dbeta = 0.0;
    const ptrdiff_t ione = 1;
    
    const double *A = mxGetPr(prhs[0]);
    const double *b = mxGetPr(prhs[1]);
    
    const ptrdiff_t mA = 3;
    const ptrdiff_t nA = 2;
    
    const ptrdiff_t mb = 2;
    const ptrdiff_t nb = 1;
    
    plhs[0] = mxCreateDoubleMatrix(mA, 1, mxREAL);
    
    double *py = mxGetPr(plhs[0]);
    
    //    N             mA   nA   alp      A  LDA  X  INCX   beta    Y   INCY
    dgemv(BLAS_NOTRANS, &mA, &nA, &dalpha, A, &mA, b, &ione, &dbeta, py, &ione);
    
    mexPrintf("%f %f %f\n", *py, *(py+1), *(py+2));
}

while simpledgemv_native use hard-coded double array in C

#include "mex.h"
#include "matrix.h"
#include "blas.h"

#if !defined(_WIN32)
#define dgemv dgemv_
#endif

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    char *BLAS_NOTRANS = "N";
    const double dalpha = 1.0, dbeta = 0.0;
    const ptrdiff_t ione = 1;
    
    const double A[] = {1.0, 0.0,
                        0.0, 1.0,
                        2.0, 0.0};
    const double *pA = A;
    const double b[] = {3.0, 4.0};
	const double *pb = b;
    
    const ptrdiff_t mA = 3;
    const ptrdiff_t nA = 2;

    const ptrdiff_t mb = 2;
    const ptrdiff_t nb = 1;
    
    plhs[0] = mxCreateDoubleMatrix(mA, 1, mxREAL);
    double *py = mxGetPr(plhs[0]);
    
    //    N             mA   nA   alp      A   LDA  X   INCX   beta    Y   INCY
    dgemv(BLAS_NOTRANS, &mA, &nA, &dalpha, pA, &nA, pb, &ione, &dbeta, py, &ione);
    
    mexPrintf("%f %f %f\n", *py, *(py+1), *(py+2));
}

So, can libmwblas be applied to do calculation on permitive double[] in C? Or I made some mistakes in writing simpledgemv_native?

Answer 1

[James Tursa]

Two problems:

2D matrices are stored column-wise by MATLAB and is assumed by the BLAS and LAPACK routines also. So this:

     const double A[] = {1.0, 0.0,
                         0.0, 1.0,
                         2.0, 0.0};

should be this instead:

    const double A[] = {1.0, 0.0, 2.0,
                        0.0, 1.0, 0.0};

If you did really want to use that first code, then you would need to reverse your A dimensions in your dgemv call (make it a 2x3) and also change your "N" to "T" in the first argument.

Then, for some reason you inserted a typo in your dgemv call, changing that last &mA to &nA (maybe you thought this was how to tell it to transpose the elements? ... it doesn't). So this

dgemv(BLAS_NOTRANS, &mA, &nA, &dalpha, pA, &nA, pb, &ione, &dbeta, py, &ione);

should be this instead:

dgemv(BLAS_NOTRANS, &mA, &nA, &dalpha, pA, &mA, pb, &ione, &dbeta, py, &ione);

[sjhstone]

So in

fullfile(matlabroot, 'extern', 'include', 'blas.h')

The definition is

/* Source: dgemv.f */
#define dgemv FORTRAN_WRAPPER(dgemv)
extern void dgemv(
    const char   *trans,
    const ptrdiff_t *m,
    const ptrdiff_t *n,
    const double *alpha,
    const double *a,
    const ptrdiff_t *lda,
    const double *x,
    const ptrdiff_t *incx,
    const double *beta,
    double *y,
    const ptrdiff_t *incy
);

I was looking into documents of Netlib's reference implementation of BLAS, in which lda is defined as "leading dimension of A".

row-major: A(i,j) = A[j + i*lda]

col-major: A(i,j) = A[i + j*lda]

My problem is perfectly solved thanks to your answer and an official support from Mathworks, but I am still curious about why "it doesn't" ...

[James Tursa]

The A matirx is always to be interpreted as being stored in column-major order. That is how MATLAB does it and that is how the BLAS and LAPACK libraries do it. If you want to use the A matrix as its transpose, you use the trans input as "T". You never change the "leading dimension of A" to accomplish this ... that is simply telling the BLAS and LAPACK libraries how to resolve memory locations for the column values.

E.g., Suppose you had a big matrix X that was 10x20 matrix stored in column-major order but you want to use a 3x5 sub-section starting at X(2,4) for the operation. In that case, you would pass 3 and 5 as the dimensions of the A matrix, you would pass the address of the X(2,4) element as the A matrix pointer (i.e., you are pointing the BLAS/LAPACK routine to start at that element), and you would pass 10 as the leading dimension of A (needed by the BLAS/LAPACK to resolve where the columns of your 3x5 submatrix are actually located in memory). If you wanted to actually multiply by the transpose of that 3x5 subsection, you would use "T" instead of "N", but you would still use 10 as the "leading dimension of A" because that is physically the number of elements between columns and the BLAS/LAPACK routines need to know this.

A consequence of this is that if you have a C/C++ native 2D array that is stored in row-major order, say 7x8, you typically need to tell the BLAS/LAPACK routines that it is a 8x7 matrix with leading dimension 8 (because that is how it is stored in memory) and use "T" for the trans to have the BLAS/LAPACK routines use it as a 7x8. (Internally the BLAS/LAPACK routines adjust memory access indexing because of the "T" to do the calculation without physically doing the transpose itself first).

[sjhstone]

This explanation is brilliant, thank you so much!