MATLAB Answers

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

2 views (last 30 days)
sjhstone
sjhstone on 9 Oct 2019
Edited: James Tursa on 15 Oct 2019
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?

Accepted Answer

James Tursa
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);
  3 Comments

Sign in to comment.

More Answers (0)

Tags

Products


Release

R2019b

Community Treasure Hunt

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

Start Hunting!