MATLAB: Use of int vs size_t in mex compilation of C-function dgemm.

blasdgemmlapackMATLABmex

I have a bit a problem understanding some issues I obtain when trying to compile a piece of C-code. In particular, with calling the dgemm function for matrix transposition.
For a bit of context: I am completely new at C-code, the code I am trying to compile is a decade old, and I had several issues compiling when running MATLAB on Windows, so I moved to Ubuntu.
I downloaded the code C-code from here, and I followed the instructions therein which in summary say:
  • Run "mex -setup" and enter the number corresponding to the template option gccopts.sh.
  • Compile with "mex fmpc_sim.c -lblas -llapack".
  • Compile with "mex fmpc_step.c -lblas -llapack".
  • Test by running the example "masses_example.
Everything compiled properly in Ubuntu, once I was able to install a version of gcc that is compatible with MATLAB, however I was never offered any options when running "mex-setup. Furthermore, although everything compiles without errors, MATLAB crashes when running the example.
Further investigation showed that crashing occurs when the code first encounters the line:
F77_CALL(dgemm)("t","n",&n,&n,&n,&fone,A,&n,eyen,&n,&fzero,At,&n);
I stripped most of the code to ensure that I can ask a more specific question. Let's say that the original C-code is called "stripped.c" and reads:
#include "blas.h"
#include "lapack.h"
#include "mex.h"
const double fone = 1;
const double fzero = 0;
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    /* problem setup */
    int i, n;
    double *dptr;
    double *A, *At;
    double *eyen;
    
    /* inputs */
    A = mxGetPr(mxGetField(prhs[0],0,"A"));
    n = (int)mxGetScalar(mxGetField(prhs[0],0,"n"));
   
    /* outputs */
    plhs[0] = mxCreateDoubleMatrix(n,n,mxREAL);
    At = mxGetPr(plhs[0]);
    eyen = malloc(sizeof(double)*n*n);
    
    /* eyen*/
    dptr = eyen;
    for (i = 0; i < n*n; i++)
    {
        *dptr = 0;
        dptr++;
    }
    dptr = dptr-n*n;
    for (i = 0; i < n; i++)
    {
        *dptr = 1;
        dptr = dptr+n+1;
    }
    /* At*/
    F77_CALL(dgemm)("t","n",&n,&n,&n,&fone,A,&n,eyen,&n,&fzero,At,&n);
    
    return;
}
Then when I execute this MATLAB-code
clear all
clc
mex -setup
mex stripped.c -lblas -llapack
sys.n=2;
sys.A=randn(sys.n);
[At]=fmpc_step_stripped(sys);
I get no compilation errors, but a matrix At filled with zeros.
If in turn I change the type declaration of variable n to:
    size_t n;
Then I get one compilation warning per each "integer" input of the F77_CALL(dgemm) function saying: "expected 'const int *' but argument is of type 'size_t *'", however At is not anymore filled with zeros, but the transpose of A as requested.
I understand the warnings, since the documentation of dgemm here says that I such inputs must be integers, but then again, it works.
Now, this is a just a very minimal working example, and many more calls to C-functions such as dgemm take place. The C-code I downloaded corresponds to what is reported in a very well cited paper, hence the compilation-running issues I am experiencing must be on my side, but I am at a lost on what could I be doing wrong. Since the C-code is so old I suspect it migth be some sort of compatibility issue, but I can't spot it also.
Any help is appreciated.

Best Answer

If you are linking to the MATLAB supplied BLAS/LAPACK libraries, then all of the integer types being used for arguments to these library functions should have the type mwSignedIndex. That is a MATLAB macro (or typedef?) that is designed to make sure your code is robust across various MATLAB versions. So if you are linking with the MATLAB supplied BLAS/LAPACK libraries, start by replacing this:
    int i, n;
        :
    n = (int)mxGetScalar(mxGetField(prhs[0],0,"n"));
with this:
    mwSignedIndex i, n;
        :
    n = (mwSignedIndex) mxGetScalar(mxGetField(prhs[0],0,"n"));
That should ensure that you are using the correctly sized integers for the BLAS/LAPACK library calls.
Also, your sample code has a memory leak for this memory:
    eyen = malloc(sizeof(double)*n*n);
I'm assuming this is because of your stripped down example and that this is not an issue in the real code? Regardless, you might want to switch from malloc/free to mxMalloc/mxFree instead.
What version of MATLAB are you using? And from the nature of your errors I am assuming 64-bit?
Related Question