C-MEX codes face problem "MATLAB has attempted to use more stack space than is available"

I used mex to convert one stiff.cpp (for example) file to stiff.mexw64 file. When I want to use this mex function stiff() in matlab, it says: Severe: MATLAB has attempted to use more stack space than is available. If a mex file was in use check for large local variables or infinite recursion.Unexpected MATLAB exception...

the code

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

void stiff(double *X,double *I,double *J,double *K)
{
int c=20;
double La=0.06;
double EA[30]={2e7,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,2e7};
double EI[30]={80,1.35,1.35,1.35,1.35,1.35,1.35,1.35,1.35,1.35,1.35,1.35,1.35,1.35,1.35,1.35,1.35,1.35,1.35,1.35,1.35,1.35,1.35,1.35,1.35,1.35,1.35,1.35,1.35,80};
double Ka[2520][2520]; 
double Kb[2520][2520]; 
double kka[84][84]; 
double kkb[84][84]; 
double ka[84][84]; 
double kb[84][84]; 
double XX[84];
int i,j,k,m,n,q,u0,ue,u1,u2,u3,u4,u5,u6,u7,u8;
double ax,ay,bx,by,ex,ey,fx,fy,a,b,e,f,dx,dy,d,A,B,C,D,E,F,AA,BB,CC,DD,EE,FF;
double k1[8][8],eye[2][2];
eye[0][0]=1;
eye[0][1]=0;
eye[1][0]=0;
eye[1][1]=1;
double  S1[8][8]={{12,0,0.54,0,-12,0,0.54,0},{0,12,0,0.54,0,-12,0,0.54},{0.54,0,0.0324,0,-0.54,0,0.0162,0},{0,0.54,0,0.0324,0,-0.54,0,0.0162}, {-12,0,-0.54,0,12,0,-0.54,0},{0,-12,0,-0.54,0,12,0,-0.54},{0.54,0,0.0162,0,-0.54,0,0.0324,0},{0,0.54,0,0.0162,0,-0.54,0,0.0324}};
double  S2[8][8]={{12,0,0.36,0,-12,0,0.36,0},{0,12,0,0.36,0,-12,0,0.36},{0.36,0,0.0144,0,-0.36,0,0.0072,0}, {0,0.36,0,0.0144,0,-0.36,0,0.0072}, {-12,0,-0.36,0,12,0,-0.36,0},{0,-12,0,-0.36,0,12,0,-0.36},{0.36,0,0.0072,0,-0.36,0,0.0144,0},{0,0.36,0,0.0072,0,-0.36,0,0.0144}};  
 for(j=0;j<30;j++)
{
  u0=0+j*4*(c+1);
  ue=4*(c+1)-1+j*4*(c+1);
  m=0;    
          for(i=u0;i<ue;i++) 
              {
                  XX[m]=X[i];
                  m++;
              }   
    for(k=0;k<20;k++)
    {
            for(i=0;i<84;i++)
        {
      for(m=0;m<84;m++)
          {
            ka[i][m]=0;
            kb[i][m]=0;
          }
        }
                u1=0+k*4;
      u2=1+k*4;
      u3=2+k*4;
      u4=3+k*4;
      u5=4+k*4;
      u6=5+k*4;
      u7=6+k*4;
      u8=7+k*4;
      ax=La*XX[u3];
                ay=La*XX[u4];
                  bx=La*XX[u7];
                by=La*XX[u8];            
                  ex=1.5*ax;
                ey=1.5*ay;
                  fx=1.5*bx;
                fy=1.5*by;
      a=sqrt(ax*ax+ay*ay);
          b=sqrt(bx*bx+by*by);
          e=sqrt(ex*ex+ey*ey);
            f=sqrt(fx*fx+fy*fy);
          dx=XX[u5]-XX[u1];
            dy=XX[u6]-XX[u2];
            d=sqrt(dx*dx+dy*dy);
  A=3/70/(La*La)*(a*a+b*b-14*La*La-6*ax*dx-6*bx*dx-6*ay*dy-6*by*dy+24*d*d);
      B=1/280/La*(b*b-a*a+2*ax*bx+2*ay*by-14*La*La-24*ax*dx-24*ay*dy+36*d*d);
  C=1/280/La*(a*a-b*b+2*ax*bx+2*ay*by-14*La*La-24*bx*dx-24*by*dy+36*d*d);
      D=1/420*(12*a*a+b*b-3*ax*bx-3*ay*by-28*La*La+3*ax*dx-3*bx*dx+3*ay*dy-3*by*dy+18*d*d);
  E=-1/840*(3*a*a+3*b*b-4*ax*bx-4*ay*by-14*La*La+6*ax*dx+6*bx*dx+6*ay*dy+6*by*dy);
      F=1/420*(a*a+12*b*b-3*ax*bx-3*ay*by-28*La*La-3*ax*dx+3*bx*dx-3*ay*dy+3*by*dy+18*d*d);
      AA=3/70/(1.5*La*1.5*La)*(a*a+b*b-14*(1.5*La*1.5*La)-6*ax*dx-6*bx*dx-6*ay*dy-6*by*dy+24*d*d);
      BB=1/280/(1.5*La)*(b*b-a*a+2*ax*bx+2*ay*by-14*(1.5*La*1.5*La)-24*ax*dx-24*ay*dy+36*d*d);
      CC=1/280/(1.5*La)*(a*a-b*b+2*ax*bx+2*ay*by-14*(1.5*La*1.5*La)-24*bx*dx-24*by*dy+36*d*d);
      DD=1/420*(12*a*a+b*b-3*ax*bx-3*ay*by-28*(1.5*La*1.5*La)+3*ax*dx-3*bx*dx+3*ay*dy-3*by*dy+18*d*d);
      EE=-1/840*(3*a*a+3*b*b-4*ax*bx-4*ay*by-14*(1.5*La*1.5*La)+6*ax*dx+6*bx*dx+6*ay*dy+6*by*dy);
      FF=1/420*(a*a+12*b*b-3*ax*bx-3*ay*by-28*(1.5*La*1.5*La)-3*ax*dx+3*bx*dx-3*ay*dy+3*by*dy+18*d*d);

        if(j==30)
        {
    k1[0][0]=AA*eye[0][0];
    k1[0][1]=AA*eye[0][1];
    k1[1][0]=AA*eye[1][0];
    k1[1][1]=AA*eye[1][1];

                k1[0][2]=BB*eye[0][0];
    k1[0][3]=BB*eye[0][1];
    k1[1][2]=BB*eye[1][0];
    k1[1][3]=BB*eye[1][1];

                k1[0][4]=-AA*eye[0][0];
    k1[0][5]=-AA*eye[0][1];
    k1[1][4]=-AA*eye[1][0];
    k1[1][5]=-AA*eye[1][1];

          k1[0][6]=CC*eye[0][0];
    k1[0][7]=CC*eye[0][1];
    k1[1][6]=CC*eye[1][0];
    k1[1][7]=CC*eye[1][1];

          k1[2][0]=BB*eye[0][0];
    k1[2][1]=BB*eye[0][1];
    k1[3][0]=BB*eye[1][0];
    k1[3][1]=BB*eye[1][1];

          k1[2][2]=DD*eye[0][0];
    k1[2][3]=DD*eye[0][1];
    k1[3][2]=DD*eye[1][0];
    k1[3][3]=DD*eye[1][1];

          k1[2][4]=-BB*eye[0][0];
    k1[2][5]=-BB*eye[0][1];
    k1[3][4]=-BB*eye[1][0];
    k1[3][5]=-BB*eye[1][1];

          k1[2][6]=EE*eye[0][0];
    k1[2][7]=EE*eye[0][1];
    k1[3][6]=EE*eye[1][0];
    k1[3][7]=EE*eye[1][1];

                k1[4][0]=-AA*eye[0][0];
    k1[4][1]=-AA*eye[0][1];
    k1[5][0]=-AA*eye[1][0];
    k1[5][1]=-AA*eye[1][1];

          k1[4][2]=-BB*eye[0][0];
    k1[4][3]=-BB*eye[0][1];
    k1[5][2]=-BB*eye[1][0];
    k1[5][3]=-BB*eye[1][1];

          k1[4][4]=AA*eye[0][0];
    k1[4][5]=AA*eye[0][1];
    k1[5][4]=AA*eye[1][0];
    k1[5][5]=AA*eye[1][1];

          k1[4][6]=-CC*eye[0][0];
    k1[4][7]=-CC*eye[0][1];
    k1[5][6]=-CC*eye[1][0];
    k1[5][7]=-CC*eye[1][1];

          k1[6][0]=CC*eye[0][0];
    k1[6][1]=CC*eye[0][1];
    k1[7][0]=CC*eye[1][0];
    k1[7][1]=CC*eye[1][1];

          k1[6][2]=EE*eye[0][0];
    k1[6][3]=EE*eye[0][1];
    k1[7][2]=EE*eye[1][0];
    k1[7][3]=EE*eye[1][1];

          k1[6][4]=-CC*eye[0][0];
    k1[6][5]=-CC*eye[0][1];
    k1[7][4]=-CC*eye[1][0];
    k1[7][5]=-CC*eye[1][1];

          k1[6][6]=FF*eye[0][0];
    k1[6][7]=FF*eye[0][1];
    k1[7][6]=FF*eye[1][0];
    k1[7][7]=FF*eye[1][1];
        n=0;
        for(i=u1;i<u8;i++)
      {
            q=0;
        for(m=u1;m<u8;m++)
        {          
          ka[i][m]=EA[j]/(1.5*La)*k1[n][q];
          kb[i][m]=EI[j]/(1.5*La*1.5*La*1.5*La)*S1[n][q];
                        q++;
        }
            n++;
      }
        }
        else
        {
          k1[0][0]=A*eye[0][0];
    k1[0][1]=A*eye[0][1];
    k1[1][0]=A*eye[1][0];
    k1[1][1]=A*eye[1][1];

          k1[0][2]=B*eye[0][0];
    k1[0][3]=B*eye[0][1];
    k1[1][2]=B*eye[1][0];
    k1[1][3]=B*eye[1][1];

          k1[0][4]=-A*eye[0][0];
    k1[0][5]=-A*eye[0][1];
    k1[1][4]=-A*eye[1][0];
    k1[1][5]=-A*eye[1][1];

          k1[0][6]=C*eye[0][0];
    k1[0][7]=C*eye[0][1];
    k1[1][6]=C*eye[1][0];
    k1[1][7]=C*eye[1][1];

          k1[2][0]=B*eye[0][0];
    k1[2][1]=B*eye[0][1];
    k1[3][0]=B*eye[1][0];
    k1[3][1]=B*eye[1][1];

          k1[2][2]=D*eye[0][0];
    k1[2][3]=D*eye[0][1];
    k1[3][2]=D*eye[1][0];
    k1[3][3]=D*eye[1][1];

          k1[2][4]=-B*eye[0][0];
    k1[2][5]=-B*eye[0][1];
    k1[3][4]=-B*eye[1][0];
    k1[3][5]=-B*eye[1][1];

          k1[2][6]=E*eye[0][0];
    k1[2][7]=E*eye[0][1];
    k1[3][6]=E*eye[1][0];
    k1[3][7]=E*eye[1][1];

          k1[4][0]=-A*eye[0][0];
    k1[4][1]=-A*eye[0][1];
    k1[5][0]=-A*eye[1][0];
    k1[5][1]=-A*eye[1][1];

          k1[4][2]=-B*eye[0][0];
    k1[4][3]=-B*eye[0][1];
    k1[5][2]=-B*eye[1][0];
    k1[5][3]=-B*eye[1][1];

          k1[4][4]=A*eye[0][0];
    k1[4][5]=A*eye[0][1];
    k1[5][4]=A*eye[1][0];
    k1[5][5]=A*eye[1][1];

          k1[4][6]=-C*eye[0][0];
    k1[4][7]=-C*eye[0][1];
    k1[5][6]=-C*eye[1][0];
    k1[5][7]=-C*eye[1][1];

          k1[6][0]=C*eye[0][0];
    k1[6][1]=C*eye[0][1];
    k1[7][0]=C*eye[1][0];
    k1[7][1]=C*eye[1][1];

          k1[6][2]=E*eye[0][0];
    k1[6][3]=E*eye[0][1];
    k1[7][2]=E*eye[1][0];
    k1[7][3]=E*eye[1][1];

          k1[6][4]=-C*eye[0][0];
    k1[6][5]=-C*eye[0][1];
    k1[7][4]=-C*eye[1][0];
    k1[7][5]=-C*eye[1][1];

          k1[6][6]=F*eye[0][0];
    k1[6][7]=F*eye[0][1];
    k1[7][6]=F*eye[1][0];
    k1[7][7]=F*eye[1][1];
        n=0;
        for(i=u1;i<u8;i++)
      {
                q=0;
        for(m=u1;m<u8;m++)
        {          
          ka[i][m]=EA[j]/La*k1[n][q];
          kb[i][m]=EI[j]/(La*La*La)*S2[n][q];
                        q++;
        }
            n++;
      }
        }    
    if(k==1)
    {
      for(i=0;i<84;i++)
      {    
        for(m=0;m<84;m++)
        {
          kka[i][m]=ka[i][m];
          kkb[i][m]=kb[i][m];
        }
      }
    }
    else
    {
      for(i=0;i<84;i++)
      {    
        for(m=0;m<84;m++)
        {
          kka[i][m]=ka[i][m]+kka[i][m];
          kkb[i][m]=kb[i][m]+kkb[i][m];
        }
      }
    }
    }

    n=0;
    for(i=u0;i<ue;i++)
      {
        q=0;
        for(m=u0;m<ue;m++)
        {          
          Ka[m][i]=kka[q][n];
          Kb[m][i]=kkb[q][n];
          q++;
        }
        n++;
      }
   }
  for(i=0;i<2520;i++)
  {
      for(j=0;j<2520;j++)
      {
          I[i*2520+j]=Ka[i][j]+Kb[i][j];
          J[i*2520+j]=Ka[i][j]+Kb[i][j];
          K[i*2520+j]=Ka[i][j]+Kb[i][j];
      }
  }
}

void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
double *mX;
double *mI;
double *mJ;
double *mK;

    plhs[0] = mxCreateDoubleMatrix(6350400,1, mxREAL);
    plhs[1] = mxCreateDoubleMatrix(6350400,1, mxREAL);
    plhs[2] = mxCreateDoubleMatrix(6350400,1, mxREAL);
    mI= mxGetPr(plhs[0]);
    mJ= mxGetPr(plhs[1]);
    mK= mxGetPr(plhs[2]);
    mX= mxGetPr(prhs[0]);
    stiff(mX,mI,mJ,mK);
}

In the C++ code, the input is a one-dimensional array with size of 2520. And the outputs include three-one-dimensional array which each size is 6350400.Is this problem caused by the big size of input and output? or others? Could you please give me some advises? THX a lot!

I've update the code. Are Ka and Kb the large local variables?

[Titus Edelhofer]

Hi,

hmm, those variables are not large. We are talking of 50 MByte variables, that shouldn't be a problem ... What do you do inside stiff, is it recursive?

Titus

Answer 1

[James Tursa]

The memory for plhs[0], plhs[1], and plhs[2] should be coming off of the heap, not the stack. So those variables should not be causing any stack problems. What else is in your code that you are not showing us (e.g., large local variables)?

[Shen Jian]

THX James! I've update the C++ code. Could you please take a look? Is there any large local variables?

[James Tursa]

Yes. These variables:

double Ka[2520][2520]; 
double Kb[2520][2520];

Each one takes over 48MB, so you've got over 96MB of stack space used by these local variables. Rewrite your code using pointers instead for these variables and use mxMalloc to allocate the memory for them (which will get the memory from the heap). E.g., to keep the same Ka[ ][ ] syntax downstream:

double *Ka[2520]; // A local array of double pointers
Ka[0] = (double *) mxMalloc(2520*2520*sizeof(double)); // Allocate one big block
for( i=1; i<2520; i++ ) {
    Ka[i] = Ka[i-1] + 2520; // Set all the other pointers
}
  :
// Then use Ka[ ][ ] syntax downstream in your code
  :
mxFree(Ka[0]); // Free the memory when done with Ka

[Shen Jian]

THX very much! James. My first mex is successful! Thank you! The code I modified:

double **Ka; 
double **Kb; 
 Ka=(double **)malloc(2520*sizeof(double*));
 Kb=(double **)malloc(2520*sizeof(double*));
for(i=0;i<2520;i++)
 {
     Ka[i]=(double *)malloc(2520*sizeof(double));
     Kb[i]=(double *)malloc(2520*sizeof(double));
 }
....../*the process of using Ka,Kb*/
 for (i=0;i<2520;i++) 
 {
        free(Ka[i]);
        free(Kb[i]);
    }
    free(Ka);
    free(Kb);

[James Tursa]

Your way also works, but my way has these advantages:

- Only one allocation call and one deallocation call (less burden on Memory Manager)
- All allocated memory for a variable is contiguous

The last one can be important if you wanted to pass the variable to a routine that expects the variable memory to be contiguous (e.g., a memcopy or similar).

Also note that by using malloc instead of mxMalloc, the burden of checking the returned pointers for NULL falls to you to avoid a potential crash. By using mxMalloc you can avoid the need for that checking.

[Shen Jian]

Dear James,

THX a lot for helping last time! The c-mex code goes well now. But, the return value Ka and Kb to MATLAB takes up a lot memory. So, I want to use Sparse matrix to return Ka and Kb. Here's the code. But it seems some errors existing again.

void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
double *mX;
mwIndex *ira,*jca,*irb,*jcb;
double *mKa;
double *mKb;
double *mS1;
double *mS2;

    plhs[0] = mxCreateSparse(u9,1,u9, mxREAL);
    plhs[1] = mxCreateSparse(u9,1,u9, mxREAL);
    mKa= mxGetPr(plhs[0]);
    ira= mxGetIr(plhs[0]);
    jca= mxGetJc(plhs[0]);
    mKb= mxGetPr(plhs[1]);
    irb= mxGetIr(plhs[1]);
    jcb= mxGetJc(plhs[1]);

    mX=mxGetPr(prhs[0]);
    mS1=mxGetPr(prhs[1]);
    mS2=mxGetPr(prhs[2]);
    stiff(mX,mS1,mS2,mKa,mKb);
}

And the value of u9 is 5080. The MATLAB help file doesn't support me. Could you please give me some advice? THX!

[James Tursa]

Sparse matrices are stored in a completely different manner than full matrices. You would have to completely re-write nearly every line of your code. Plus, inserting a value at an arbitrary location is not trivial. If you really want to explore this option let me know and I can work up an example, but be forewarned it will be a LOT of work for you.

[Shen Jian]

Thanks for replying! I think I can find other ways to accelerate my code in MATLAB. THX again!