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

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.

Answer 1

[James Tursa]

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?

[Bernardo Hernandez]

Dear James,

Thanks for your answer.

I am using R2017b 64-bit. I also assume I am using the MATLAB supplied BLAS and LAPACK libraries because I haven't done anything to actively replace them.

When replacing "int" by "mwSignedIndex" I see a similar behaviour as to when I used "size_t". That is, when compiling I get a warning for every integer entry on the dgemm call, saying: "expected 'const int *' but argument is of type 'mwSignedIndex *'". Nevertheless, the compilation is succesful and the code runs as expected.

I know that I only posted a stripped down version of the code, but could this be an indicator that the sole reason the C-code I downloaded doesn't compile is version compatibility issues?

With respect ot the "malloc" memory leak. No, this is not my doing and it is all over the code. I will try to replace it with what you suggest, provided I don't break anythinng else in the process.

I will see if I can extrapolate your answer to get the entire code running. If not, I will post another question, as I assume this one should be considered answered.

Thank you very much.

[James Tursa]

"... "expected 'const int *' but argument is of type 'mwSignedIndex *'". Nevertheless, the compilation is succesful and the code runs as expected ..."

This should not be happening on a 64-bit system. That is, the libraries definitely use 64-bit integers (mwSignedIndex) but your compiler thinks they use 32-bit integers (int). I can only postulate that there is a bad prototype involved somewhere incorrectly telling the compiler that those arguments are int pointers. So, maybe you are including the wrong blas.h file that is from a 32-bit install? Is the blas.h file you are using actually included with your package? If so, don't use it ... use the one from your 64-bit install instead. Bottom line is to make sure you are using the blas.h file that comes with your MATLAB version, as it should have the correct prototypes for the BLAS function arguments. (You could also check that the blas.h file for your MATLAB version has correct prototypes, but I think it unlikely that this is the problem.)

"... I know that I only posted a stripped down version of the code, but could this be an indicator that the sole reason the C-code I downloaded doesn't compile is version compatibility issues? ..."

Yes.

"... With respect ot the "malloc" memory leak. No, this is not my doing and it is all over the code. I will try to replace it with what you suggest ..."

The main problem is that there is no free(eyen) call at the end of the routine to free the memory allocated by malloc( ), so there will be a permanent memory leak for this that will not be recovered until you restart MATLAB. Every malloc( ) call needs to be paired with a free( ) call. By using mxMalloc( ) instead of malloc( ), you at least get some protection against permanent memory leaks with help from the MATLAB Memory Manager if you don't have corresponding mxFree( ) calls. But really, for good programming practice, you should have those mxFree( ) calls in the code.

Answer 2

[Bernardo Hernandez]

Dear James,

First let me tell you that I was able to run the code following your first set of suggestions, so thank you very much.

With respect to your additional comments:

1. I was indeed using libraries provided with the C-code I downloaded. I removed them and copied in the blas.h and lapack.h files from the MATLAB installation folder. This, however, resulted in the stripped version of the code not compiling (it just stays in busy for too long). I am not sure what could the problem be, but in any case I am satisfied with the fact that the code runs.
2. Thanks for the memory leak suggestion, I will make sure to include it in the code to make it more robust.

Kind regards.

[James Tursa]

"... I was able to run the code following your first set of suggestions ..."

That's good, but the danger of using incorrect prototypes in general is that an incorrect automatic type conversion for an argument could be made by the compiler. This won't happen in your particular case since all of your argument issues are with pointers (pointer-to-int is the same size as pointer-to-mwSignedIndex), which is why the code still works with an incorrect prototype.

"... This, however, resulted in the stripped version of the code not compiling (it just stays in busy for too long) ..."

This sounds like a new problem ... why wouldn't the compile at least finish? Are your include files recursive?

[rupal]

hi i am also using this code using all your suggesstions I am able to complied the c file but when i am trying to run the examples still my matlab crashed. If your code is runnig (as well as example) than please share with me. This will be very helpfull for me.