2 views (last 30 days)

Show older comments

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?

James Tursa
on 10 Oct 2019

Edited: James Tursa
on 10 Oct 2019

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

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!