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

조회 수: 1(최근 30일)
sjhstone 2019년 10월 9일
편집: James Tursa 2019년 10월 15일
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);
const double *b = mxGetPr(prhs);
const ptrdiff_t mA = 3;
const ptrdiff_t nA = 2;
const ptrdiff_t mb = 2;
const ptrdiff_t nb = 1;
plhs = mxCreateDoubleMatrix(mA, 1, mxREAL);
double *py = mxGetPr(plhs);
// 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 = mxCreateDoubleMatrix(mA, 1, mxREAL);
double *py = mxGetPr(plhs);
// 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 2019년 10월 10일
편집: James Tursa 2019년 10월 10일
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};
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);
dgemv(BLAS_NOTRANS, &mA, &nA, &dalpha, pA, &mA, pb, &ione, &dbeta, py, &ione);
##### 댓글 수: 3표시숨기기 이전 댓글 수: 2
sjhstone 2019년 10월 11일
This explanation is brilliant, thank you so much!

댓글을 달려면 로그인하십시오.

R2019b

### Community Treasure Hunt

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

Start Hunting!