MEX Crash due to dynamic memory allocation

This is my first time attempting to convert a MATLAB function into a C++ executable. I am having difficulty dynamically allocating space to arrays that are created and used only within the MEX function. This is necessary because the function will be used in an iterative process in which the number of points of interest for the calculations will change. I am able to compile the code successfully in its current state, however,MATLAB crashes without any specific error when the function is run.

The purpose of this code is to calculate the stresses [s11, s22, s12] occurring between different points within a 2D domain. There are 8 coordinate frames of reference to be considered. The distance values passed in are already resolved into the reference frame of the specific point of interest. I have attache the entire script below with comments for clarity. I have also attached a Matlab script that should generate the inputs necessary for the MEX file.

#include <math.h>
#include <matrix.h>
#include <mex.h>
#include <iostream>

# define PI 3.1415926535897932

//MATLAB function definition
//[stress_11,stress_22,stress_12] = SMex(dXi_Rx,dXi_Ry,b_vec2,slip_sel,theta,h,mu,w,l,nu)

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	//declare variables
    mxArray *g_out, *h_out, *i_out;
    const mwSize *dims;
    double *X_p, *Y_p, *b_p, *theta, *slip_p, *mu, *h_d, *nu, *w, *l, *s11, *s22, *s12;
    int dimx, dimy, *slip;
    int j,k,q,h,z, aridx, sidx;

//associate inputs
	mxArray *dupArray[10];
	for(aridx = 0; aridx < 10; ++aridx){
		dupArray[aridx] = mxDuplicateArray(prhs[aridx]);
	}

//figure out dimensions
    dims = mxGetDimensions(prhs[0]);
    dimy = (int)dims[0]; dimx = (int)dims[1];

//associate outputs
    g_out = plhs[0] = mxCreateDoubleMatrix(1,dimx,mxREAL);
    h_out = plhs[1] = mxCreateDoubleMatrix(1,dimx,mxREAL);
	i_out = plhs[2] = mxCreateDoubleMatrix(1,dimx,mxREAL);

//associate pointers
	 X_p = mxGetPr(dupArray[0]);
    Y_p = mxGetPr(dupArray[1]);
    b_p = mxGetPr(dupArray[2]);
    slip_p = mxGetPr(dupArray[3]);
	theta = mxGetPr(dupArray[4]);
    h_d = mxGetPr(dupArray[5]);
	mu = mxGetPr(dupArray[6]);
	nu = mxGetPr(dupArray[7]);
	w = mxGetPr(dupArray[8]);
	l = mxGetPr(dupArray[9]);
    s11 = mxGetPr(g_out);
    s22 = mxGetPr(h_out);
	s12 = mxGetPr(i_out);

//Convert h_d from double to int for for loop
	h = (int)h_d[0];

//Convert slip system values from double to int for switch structure
	for(z=0;z<dimy;z++){
		slip[z] = (int)slip_p[z];
	}

//Declare Variables to be used in Stress calculations
//Dynamically allocate space as the number of points of interest will change 
	double term1, num, deno, s_term[3];
	double b, X, Y;
	int Ytemp, index;

	double *s[3][8]; 
	for(sidx = 0; sidx < 24; sidx){
		s[sidx] = new double [dimx];
	}

// Stresses are calculated at each dislocation, 11 identical domains are considered stacked ontop of one another.
    for(j=0;j<=(2*h);j++)
    {
		
		Ytemp =(j-h);
		for(k=0;k<dimx;k++)
		{
			for(q=0;q<dimy;q++)
			{
				index = k*dimy+q;

				if(k == q){
					s_term[0] = 0;
					s_term[1] = 0;
					s_term[2] = 0;
				}
				else{
				b = b_p[q];
				X = X_p[index]/(w[0]/2);
				Y = Y_p[index]/(w[0]/2) + Ytemp;

				term1 = mu[0]*PI*b/(2*PI*(1-nu[0])*w[0]);
				num = 1-cos (PI*X)*cosh (PI*Y);
				deno = 1/(cosh (PI*Y) - cos (PI*X));

				s_term[0] = -1*term1*deno*(2*sinh (PI*Y) + PI*Y*num*deno);
				s_term[1] = 0;
				s_term[2] = 0;
				}
// Stresses aributted to dislocations corresponding to the same slip system are summed
                              for (sidx = 0; sidx < 3; ++sidx){
			            qidx = slip[q];
				    s[sidx][qidx][k] = s[sidx][qidx][k] + s_term[sidx];
				}
				
			}
                 }
    }
// Stresses are resolved into the reference coordinate frame and returned to MATLAB
		for(k=0;k<dimx;k++){
			for(q=0;q<8;q++){
			s11[k] = s11[k] +cos(theta[q])*cos(theta[q])*s[0][q][k] + sin(theta[q])*cos(theta[q])*2*s[2][q][k] + sin(theta[q])*sin(theta[q])*s[1][q][k];

			s22[k] = sin(theta[q])*sin(theta[q])*s[0][q][k] - cos(theta[q])*sin(theta[q])*2*s[2][q][k] + cos(theta[q])*cos(theta[q])*s[1][q][k] ;
				     
			s12[k] = -cos(theta[q])*sin(theta[q])*s[0][q][k] - (cos(theta[q])*cos(theta[q]) - sin(theta[q])*sin(theta[q]))*s[2][q][k] + cos(theta[q])*sin(theta[q])*s[1][q][k];
			}
		}
    return;

}

[Guillaume]

It's not the allocations that are going to cause a crash. it's most likely because you use something that has not been properly allocated, or one of your pointer is invalid.

You certainly like copy/pasting. Just as in matlab, numbered variables in C (and any other language for that matter) are always a bad idea. Wouldn't

mxArray* dupArray[10];
for (aridx = 0; ardix < 10; ++aridx)
    dupArray[aridx] = mxDuplicateArray(prhs[aridx]);

be a lot simpler than your bunch of xx_in?

Same for your wall of sxx_x arrays. If you'd declared

double* s[3][8];
for (scol = 0; scol < 3; ++scol)
    for (srow = 0; srow < 8; ++srow)
    s[scol, srow] = (double*) MxMalloc(dimx * sizeof(double));  %is MxMalloc necessary

then, your huge case is just:

for (sidx = 0; sidx < 3; ++sidx){}
    qidx = (int)slip[q];
    s[sidx][qidx][k] += s[sidx][qidx][k];
}

However, why are you duplicating the input arrays? As far as I can tell you don't modify the inputs.

Why are all these s arrays allocated with mxMalloc. If you were actually using C++ (your code is just C at the moment) you could use new to ensure that your arrays are deallocated at the end of the function. At the moment, I can't see them being deallocated.

More importantly, I don't see you assigning anything to the output arrays.

edited to fix several mistakes as pointed below

[James Tursa]

@Guillaume: Typo, I think you meant this:

mxArray *dupArray[10];

Also, this won't compile since type type of pointer on the lhs doesn't match the type of pointer on the rhs:

double* s[3][8];
    :
    s[sidx] = (double*) etc.

[morri295]

Thank you for your input. This is my first time coding in C/C++ so I am not comfortable with the syntax. I used your first two snippets of code as that does look much better.

I am duplicating the input arrays because that is what I saw others do when making MEX files for MATLAB. I assumed it was necessary.

I don't think your third piece of code will do what is desired. The array "s" would only have 24 elements but dimx and dimy are signifcantly larger than this.

The case structure is defined by the variable slip. each row of the 2D array is assigned a specific reference frame as defined by slip. I only want to add values in the same column that correspond to the same slip value.

[Guillaume]

@James, @morri

Well, it's been over 10 years since I wrote any C++, exact syntax is becoming hazy... Follow the spirit of the comment (also why i didn't post it as an answer), not the actual riddled with typos and oversight code.

@morri,

The idea with the 2nd and 3rd piece is to create a 3x8 array of pointers to array, so s[i][j] replaces your si_j and s[i][j][k] replaces your si_j[k]. The only thing is I forgot that the j index in your case is not q but slip[q], so:

for (sidx = 0; sidx < 3; ++sidx){}
    qidx = (int)slip[q];
    s[sidx][qidx][k] += s[sidx][qidx][k];
}

Oh, and I forgot to mention, the spacing in your

dimx* sizeof(*s11_1)  %inside the MxMalloc calls

is extremely misleading. The first * is a multiplication but with that spacing it looks like a pointer declaration. The compiler is going to understand it fine, but you're going to confuse human readers.

[morri295]

Thank you for the comment, that is a much more concise method of writing the code. I am still having issues with the code crashing when I try to execute it. It will compile but not run.

Answer 1

[Robert Lummert]

Your function never checks nrhs>=1 but accesses prhs[0]. If called with no input argument, it should crash when accessing the non-existant input in line 29 already.

[morri295]

I checked the nrhs. It gives me an integer of 10. To my understanding, I should then be able to access prhs[1] up to prhs[9].

[Robert Lummert]

Well, prhs[0] up to prhs[9], but yes, in that case the use of prhs shouldn't cause the crash.