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

blasMATLABmex

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?

Best Answer

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);
Related Question