Help to debug mexFunction

I have a very simple mexFunction but it still causes Matlab to crash. Can someone help to debug where the error is?

#include "geometryMatrix64.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
  mxArray *inputLat;
  double *arrayLat= NULL;
  double *outputGeo= NULL;
  int outputRowLen = 0, outputColLen = 0, i = 0, j = 0;

  inputLat = prhs[0];

  arrayLat = mxGetPr(inputLat);

  mexPrintf("Starting...............");
  outputRowLen = 3;
  outputColLen = 1;
  plhs[0] = mxCreateDoubleMatrix(outputRowLen, outputColLen, mxREAL);
  outputGeo = mxGetPr(plhs[0]);
  for (i=0;i<outputRowLen;i++){
    for(j=0; j<outputColLen; j++){
      outputGeo[i*outputRowLen+j] = 0;
    }
  }

  mexPrintf("Calculation of geometry matrix done!");

  
  return ;
}

I compiled with the following matlab code:

for i=1:10
    i
    clear all;
    load('parameters.mat')
    blaslib = fullfile('C:\Program Files\MATLAB\R2011b\extern\lib\win64\microsoft\libmwblas.lib');
    lapacklib = fullfile('C:\Program Files\MATLAB\R2011b\extern\lib\win64\microsoft\libmwlapack.lib');
    mex('-largeArrayDims', '-compatibleArrayDims', 'geometryMatrix64.c', 'mathFunction64.c', blaslib, lapacklib)
    %[Geoc] = geometryMatrix64(lat, long, incl, peri, seax, GM, r_E, ele_mask, time, theta);
    [Geoc] = geometryMatrix64(lat);
end;

On the command window, such message is shown and crashed:

i = 
    1
Starting...............Calculation of geometry matrix done!

i =
    2
Starting...............Calculation of geometry matrix done!

Restarted computer and I got the same result but the iteration went up to i = 3.

I'm sorry for posting so often lately about mex functions. But I'm really stuck and am getting more and more frustrated since a simple mexFunction like this can't even work.

[Jane Jean]

Just to add a piece of additional information.

I'm running Matlab R2011b 64-bit on a 64-bit Windows 7.

Answer 1

[Jan]

I do not see a reason to recompile the function repeatedly.

outputRowLen = 3;
outputColLen = 1;
plhs[0] = mxCreateDoubleMatrix(outputRowLen, outputColLen, mxREAL);

Now plhs[0] is a [3 x 1] vector.

for (i=0;i<outputRowLen;i++){
  for(j=0; j<outputColLen; j++){
    outputGeo[i*outputRowLen+j] = 0;
  }
}

But here the highest index is i=2 and j=0 => [i*outputRowLen+j] = 2*3+0 = 6. Therefore you are writing in unallocated memory. You confused Row with Col in these loops.

[Jane Jean]

Thank you very much. I was too desperate that I got blinded by my impatience. :) You saved my day!

[Jan]

You are really welcome! I've struggled too often with the [i*Len+j] bugs in my own code also. Now I'm using linear indexing always:

for (i=0; i<RowLen*ColLen; i++) {X[i] = ...}

or:

Xf = X + RowLen*ColLen; while(X < Xf) {*X++ = ...}

This does not necessarily avoid bugs or increase the readability. But at least I cannot confuse the columns and rows anymore...

[James Tursa]

@Jane: FYI, as a general rule, to make better use of cache memory it is best to drag your arrays through the cache in linear fashion. Since MATLAB stores matrices in column order, that means you should have the outer loop be the columns and the inner loop be the rows. Also note that I have changed the indexing calculation. E.g.,

for( j=0; j<outputColLen; j++ ){

for( i=0; i<outputRowLen; i++ ){

outputGeo[i + outputRowLen*j] = 0;

}

}

The way you currently have the loops ordered you are jumping around in memory as you drag the matrix through the high level cache, which degrades performance. You might not notice a difference for small sized arrays, but you will definitely notice it once you start working with large arrays. Best to get in the habit of coding your loops in this order from the beginning. This also makes it easier to parallelize your code in the future, and has the benefit that it can decrease the amount of memory that is required to be loaded into the individual cores.