MexFunction Crash when running two for-loops

Question

0 个投票

I have written a mexFunction with two for-loops. The strange thing is that, if I comment one of the loops, the mexFunction is called successfully. However, if I let both loops run, Matlab crashes. I use Matlab 2017a.

What is the problem? Memory issues?

Yutao

#if !defined(_WIN32)
#define dgemm dgemm_
#define dgemv dgemv_
#endif
#include "mex.h"
#include "blas.h"
/*
 *r: block row index
 *c: block column index
 *Gi: block data that is filled in
 *rows: number of rows of the block Gi
 *cols: number of columns of the block Gi
 *G: the original matrix
 *N: the number of row blocks in G
 *
 *This function assumes all blocks in G have the same dimension
*/
void Block_Fill(int r, int c, double *Gi, int rows, int cols, double *G, int N){
      int i,j,spo,spi;
      spo = r*rows+c*cols*N*rows;
      for(j=0;j<cols;j++){
          spi=spo+ j*N*rows;
          for(i=0;i<rows;i++){
              G[spi+i]=Gi[i+j*rows];
          }
      }
  }
void Block_Access(int r, int c, double *Gi, int rows, int cols, double *G, int N){
      int i,j,spo, spi;
      spo = r*rows+c*cols*N*rows;
      for(j=0;j<cols;j++){
          spi=spo+ j*N*rows;
          for(i=0;i<rows;i++){
              Gi[i+j*rows] = G[spi+i];
          }
      }
  }
void
mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
{    
    /*Inputs*/
    // A{} prhs[0]
    // B{} prhs[1]   
    double *ds0 = mxGetPr(prhs[2]);
    double *a = mxGetPr(prhs[3]);        
    int nx = mxGetScalar(prhs[12]);
    int nu = mxGetScalar(prhs[13]);
    int N = mxGetScalar(prhs[16]);
      int nz = nx+nu;
      /*Outputs*/
      double *G, *L;
      plhs[0] = mxCreateDoubleMatrix(N*nx, N*nu, mxREAL);
      G = mxGetPr(plhs[0]);    
      plhs[1] = mxCreateDoubleMatrix((N+1)*nx, 1, mxREAL);
      L = mxGetPr(plhs[1]);   
      /*Allocate memory*/
      mwIndex i,j;
      const mxArray *cell_element;  
      double *Ai, *Bi; //from cells
      double *Gi, *Ci, *Li, *ai; // Intermediate data
      Gi = mxGetPr(mxCreateDoubleMatrix(nx,nu,mxREAL));
      Ci = mxGetPr(mxCreateDoubleMatrix(nx,nu,mxREAL));
      Li = mxGetPr(mxCreateDoubleMatrix(nx,1,mxREAL));
      ai = mxGetPr(mxCreateDoubleMatrix(nx,1,mxREAL));
      /*Initialization*/
      memcpy(&ai[0],&ds0[0], nx*sizeof(double));  
      char *nTrans = "N", *Trans="T";
      double one = 1.0, zero = 0.0;
      int ONE = 1;
      /*Start the loop*/
      // compute G
      for(i=0;i<N;i++){
          cell_element = mxGetCell(prhs[1], i);
          Bi = mxGetPr(cell_element);
          Block_Fill(i,i,Bi,nx,nu,G,N);     
          for (j=i+1;j<N;j++){
              Block_Access(j-1,i,Gi,nx,nu,G,N);
              cell_element = mxGetCell(prhs[0], j);
              Ai = mxGetPr(cell_element);
              dgemm(nTrans, nTrans, &nx, &nu, &nx, &one, Ai, &nx, Gi, &nx, &zero, Ci, &nx);
              Block_Fill(j,i,Ci,nx,nu,G,N);
          }
      }
      // compute L
      for(i=0;i<N;i++){       
          memcpy(&Li[0],&ai[0], nx*sizeof(double));
          memcpy(&L[i*nx],&Li[0], nx*sizeof(double));
          memcpy(&ai[0],&a[i*nx], nx*sizeof(double));
          cell_element = mxGetCell(prhs[0], i);       
          Ai=mxGetPr(cell_element);      
          dgemv(nTrans,&nx,&nx,&one,Ai,&nx,Li,&ONE,&one,ai,&ONE);                    
      }
      memcpy(&L[N*nx],&ai[0], nx*sizeof(double));
      return;
}

2 个评论
显示无隐藏无

James Tursa 2017-9-7

编辑：James Tursa 2017-9-7

"... I have written a mexFunction with two for-loops. ..."

I count seven.

"... If I comment one of the loops, the mexFunction is called successfully. However, if I let both loops run, Matlab crashes. ..."

Do you think that maybe, just maybe, your readers might be interested to know which loop seems to be the problem?

Also, don't use int for the integer arguments to the TMW supplied BLAS or LAPACK functions. Always use mwSignedIndex instead.

Yutao Chen 2017-9-22

Thank you James,

When I say two for-loops I mean two outer loops.

I did follow your suggestion to change each int to mwSize or mwIndex. The code works! It has bothered me a month, but I don't expect the error is like this.

Thank you again

Yutao

请先登录，再进行评论。

请先登录，再回答此问题。

Follow Question