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

1 visualizzazione (ultimi 30 giorni)
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?

Risposta accettata

James Tursa
James Tursa il 10 Ott 2019
Modificato: James Tursa il 10 Ott 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 Commenti
James Tursa
James Tursa il 10 Ott 2019
Modificato: James Tursa il 15 Ott 2019
The A matirx is always to be interpreted as being stored in column-major order. That is how MATLAB does it and that is how the BLAS and LAPACK libraries do it. If you want to use the A matrix as its transpose, you use the trans input as "T". You never change the "leading dimension of A" to accomplish this ... that is simply telling the BLAS and LAPACK libraries how to resolve memory locations for the column values.
E.g., Suppose you had a big matrix X that was 10x20 matrix stored in column-major order but you want to use a 3x5 sub-section starting at X(2,4) for the operation. In that case, you would pass 3 and 5 as the dimensions of the A matrix, you would pass the address of the X(2,4) element as the A matrix pointer (i.e., you are pointing the BLAS/LAPACK routine to start at that element), and you would pass 10 as the leading dimension of A (needed by the BLAS/LAPACK to resolve where the columns of your 3x5 submatrix are actually located in memory). If you wanted to actually multiply by the transpose of that 3x5 subsection, you would use "T" instead of "N", but you would still use 10 as the "leading dimension of A" because that is physically the number of elements between columns and the BLAS/LAPACK routines need to know this.
A consequence of this is that if you have a C/C++ native 2D array that is stored in row-major order, say 7x8, you typically need to tell the BLAS/LAPACK routines that it is a 8x7 matrix with leading dimension 8 (because that is how it is stored in memory) and use "T" for the trans to have the BLAS/LAPACK routines use it as a 7x8. (Internally the BLAS/LAPACK routines adjust memory access indexing because of the "T" to do the calculation without physically doing the transpose itself first).
sjhstone
sjhstone il 11 Ott 2019
This explanation is brilliant, thank you so much!

Accedi per commentare.

Più risposte (0)

Tag

Prodotti


Release

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by