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

조회 수: 1(최근 30일)
sjhstone
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[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
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};
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
sjhstone
sjhstone 2019년 10월 11일
This explanation is brilliant, thank you so much!

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

추가 답변(0개)

태그

제품


릴리스

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by