Convert Matlab function to Mex file - For loop conversion

I am trying to optimize my code by converting the following for loop to mex file. However, I am having a few problems converting them. Can anyone help me ?

    Ixyz1 =  zeros(xgridsize * ygridsize * zgridsize, 1,'single');
    Ixyz2 =  zeros(xgridsize * ygridsize * zgridsize, 1,'single');
    Ixyz3 =  zeros(xgridsize * ygridsize * zgridsize, 1,'single');
    Ixyz4 =  zeros(xgridsize * ygridsize * zgridsize, 1,'single');
    Ixyz5 =  zeros(xgridsize * ygridsize * zgridsize, 1,'single');
    Ixyz6 =  zeros(xgridsize * ygridsize * zgridsize, 1,'single');
    Ixyz7 =  zeros(xgridsize * ygridsize * zgridsize, 1,'single');
    Ixyz8 =  zeros(xgridsize * ygridsize * zgridsize, 1,'single');
    Ixyz9 =  zeros(xgridsize * ygridsize * zgridsize, 1,'single');
    Ixyz10 = zeros(xgridsize * ygridsize * zgridsize, 1,'single');

    for i = 1:1:length(Ixyz1)   
        Ixyz1(i)  = Ixy(timedifference1(i) + 2001,1);
        Ixyz2(i) = Ixy2(timedifference2(i) + 2001,1);
        Ixyz3(i) = Ixy3(timedifference3(i) + 2001,1);
        Ixyz4(i) = Ixy4(timedifference4(i) + 2001,1);
        Ixyz5(i) = Ixy5(timedifference5(i) + 2001,1);
        Ixyz6(i) = Ixy6(timedifference6(i) + 2001,1);
        Ixyz7(i) = Ixy7(timedifference7(i) + 2001,1);
        Ixyz8(i) = Ixy8(timedifference8(i) + 2001,1);
        Ixyz9(i) = Ixy9(timedifference9(i) + 2001,1);
        Ixyz10(i) = Ixy10(timedifference10(i) + 2001,1);
    end

timedifference1 to timedifference10 have the same size as Ixyz1 to Ixyz10.

xgridsize , ygridsize are usually greater than 200. and zgridsize is around 100. Ixy to Ixy10 have the same size of 4000*1.

[Walter Roberson]

Please say more about the difficulty you have in the conversion?

Is there a reason you did not use vectorized code

Ixyz1 = Ixy(timedifference1 + 2001);
Ixyz2 = Ixy2(timedifferenc2 + 2001);

and so on?

[aa39998]

I tried vectorizing the code. It didn't show much significant improvement in terms of run time. The reason being that the size of timedifference1 is too large, almost greater than 2^21. So for loop and vectorization have roughly the same performance under this condition.

As for the conversion, I created the matrix for return from Ixyz1 to Ixyz10.When I compile the mex file, it compiled successfully. However, when I run my script, matlab crashes immediately. Thus, I have no idea what happened to my conversion. One of my speculation is that the mxCreateDoubleMatrix I created is too large for the memory to hold. But, I can't get away with creating a huge matrix after my computation in C. Is there any one I can get around with it ?

This is the code I converted : (It only has one Ixyz1 return, but the idea is the same if I want all Ixyz1 matrices to be returned )

#include "mex.h"
void fillgrid(double Ixyzlength, double Ixy[], double A[], double *Ixyz1){
  for (int i = 0 ; i <Ixyzlength; i++) {
      *(Ixyz1+i)=Ixy[(int)A[i]+2000];
  }
}
void mexFunction(int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[]){
double *A, *Ixy, *Ixyz1;
double Ixyzlength;
/* Check for proper number of arguments. */
  if (nrhs != 3) 
  {
     mexErrMsgIdAndTxt("MATLAB:fillgrid:invalidNumInputs","Incorrect number of inputs");
  }
  if (nlhs > 2) 
  {
     mexErrMsgIdAndTxt("MATLAB:fillgrid:maxlhs", "Too many output  arguments.");
  }
/* Create matrix for the return argument. */
  plhs[0] = mxCreateDoubleMatrix((int) Ixyzlength, 1,mxREAL);
/* Assign pointers to each input and output. */
  A = mxGetPr(prhs[0]);
  Ixy = mxGetPr(prhs[1]);
  Ixyzlength = mxGetScalar(prhs[2]);
//Output
  Ixyz1 = mxGetPr(plhs[0]);
/* Call the subroutine. */
  fillgrid(Ixyzlength, Ixy , A, Ixyz1);
  mxDestroyArray(Ixyz1);
}

Answer 1

[OCDER]

1st issue is that you called mxCreateDoubleMatrix BEFORE Ixyzlength was initialized. Also, mxCreateDoubleMatrix 1st input takes type mwSize (which is of size_t, not int), though it'll cast it. Careful of the casting rules.

 //ERROR. Ixyzlength is declared but not initialized.
  double Ixyzlength;  
  plhs[0] = mxCreateDoubleMatrix((int) Ixyzlength, 1,mxREAL);  

 //FIX by moving mxCreateDoubleMatrix after you get Ixyzlength.
  double Ixyzlength = mxGetScalar(prhs[2]);
  plhs[0] = mxCreateDoubleMatrix((mwSize) Ixyzlength, 1,mxREAL);

2nd issue is that you are destroying your Ixyz1 output array, which you shouldn't as stated here: https://www.mathworks.com/help/matlab/apiref/mxdestroyarray.html

 //DELETE THIS
  mxDestroyArray(Ixyz1);

3rd issue, it seems that Ixyzlength is just the size of Ixy. You could use mxGetM , mxGetN , or mxGetNumberOfElements to get this information as type mwSize. That way, you don't have to provide a 3rd input to your mex function.

4th issue, do not label variables Ixyz1/2/3/4/5/6/7/9/10 . That'll result in a huge code that's hard to modify. Instead, use cell arrays in your case so you can do this:

%Ixy is also converted to a cell array
Ixyz = repmat({zeros(xgridsize * ygridsize * zgridsize, 1,'single')}, 10, 1);
for j = 1:length(Ixyz)
    for i = 1:length(Ixyz{1})
       Ixyz{j}(i) = Ixy{j}(timedifference1(i) + 2001,1)
    end
end

[aa39998]

Thanks for the answer. I tried fixing the issue that you mentioned. However, I still got the segmentation violation after running the mex file in my script.What could have caused the issue ?

Also, Ixyzlength is the size of A , not Ixy. Just to clarify.

[James Tursa]

Please post your current C code, and also the m-code that you are using to call the C code.

[aa39998]

The input arugments for the c code are 2 matrices. Matrix A is a int16 matrix of size 2,302,901 *1. Matrix Ixy is a type single, matrix of size 4001*1. The output should be only 1 matrix ,Ixyz , type double , size 2302901*1.

In my matlab script. I am calling this function as

fillgrid(A, Ixy)

The Current C code is as follows :

 #include "mex.h"
void fillgrid(int Ixyzlength, float Ixy[], int A[], double *Ixyz1){
  for (int i = 0 ; i <  Ixyzlength; i++) {
    *(Ixyz1+i)=Ixy[(int)A[i]+2000];
  }
}
void mexFunction(int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[]) {
  double *Ixyz1;
  int *A;
  float *Ixy;
/* Check for proper number of arguments. */
    if (nrhs != 2) {mexErrMsgIdAndTxt("MATLAB:fillgrid:invalidNumInputs", "Incorrect number of inputs");
  }
    if (nlhs > 2) {mexErrMsgIdAndTxt("MATLAB:fillgrid:maxlhs", "Too many output arguments.");
  }
    int Ixyzlength = mxGetNumberOfElements(prhs[0]);
/* Create matrix for the return argument. */
    plhs[0] = mxCreateDoubleMatrix((mwSize) Ixyzlength, 1,mxREAL);
/* Assign pointers to each input and output. */
    A = (int*) mxGetData(prhs[0]);
    Ixy = (float*) mxGetData(prhs[1]);
//Output
    Ixyz1 = mxGetPr(plhs[0]);
/* Call the subroutine. */
    fillgrid(Ixyzlength, Ixy , A, Ixyz1);
}

[James Tursa]

If "A" is an int16, then this

int *A;

should be this instead

short *A;

And of course all of the code that uses "A" needs to be adjusted as well. E.g., the (int *) cast should be (short *), the fillgrid "int A[]" needs to be "short A[]", etc.

Or, if you want to use the definitions from the tmwtypes.h file, you could use INT16_T instead of short.

You need to put in many more argument checks in your mex code to ensure that the inputs are exactly the type & size etc as expected.