MATLAB: MexGetPr error when freeing the data

Question

Hi I am using mex function to write a c code with mkl library on linux. I use the following makefile:

MATLAB = /usr/local/packages/matlabR2016a
MKL = /u/c/radhiali/intel/mkl
CFLAGS   = -Wall -g -O -shared -fPIC -fexceptions -fno-omit-frame-pointer \
           -pthread -std=gnu99
CPPFLAGS = -DMX_COMPAT_32 -DMKL_ILP64 -DMATLAB_MEX_FILE -DNDEBUG \
           -I$(MATLAB)/extern/include -I$(MATLAB)/simulink/include \
           -I$(MKL)/include
LDFLAGS  = -Wl,--version-script,$(MATLAB)/extern/lib/glnxa64/mexFunction.map \
           -L$(MATLAB)/bin/glnxa64 -lmx -lmex -lmat -lm -lstdc++ \
           -Wl,--start-group \
           $(MKL)/lib/intel64/libmkl_intel_ilp64.a \
           $(MKL)/lib/intel64/libmkl_sequential.a \
           $(MKL)/lib/intel64/libmkl_core.a \
           -Wl,--end-group
mexDAFEM.mexa64: mexDAFEM.c
        $(CC) $(CFLAGS) $(CPPFLAGS) -o $@ $^ $(LDFLAGS)

I have included part of my function below. It seems the error comes from the mkl_free statements of pointers r and rlist in the function neighbours. The error seems to be a segmentation fault. It worked for windows but not for linux. Could smeone help me to fix this problem? The variables are obtained by a mxGetPr function. Thank you and have a good one.

Regards Ali Radhi

#include <mex.h>
#include <stdio.h>
#include <math.h>
#include <matrix.h>
#include "mkl_spblas.h"
#include "mkl_types.h"
#include "mkl.h"
#include "mkl_pblas.h"
#include "mkl_pardiso.h"
#include "mexNeighbour.h"
#include "mexResidual.h"
#include "mexStiffness.h"
void pbc (double *posT, double *pcell, MKL_INT Na) 
{
    MKL_INT ii;

      for(ii=0;ii<3*Na;ii+=3){
          if (posT[ii]>pcell[0]){
              posT[ii]=posT[ii] - 2*pcell[0];
          }
          else if (posT[ii]<-pcell[0]){
              posT[ii]=posT[ii] + 2*pcell[0];
          }
      }

      for(ii=1;ii<3*Na;ii+=3){
          if (posT[ii]>pcell[1]){
              posT[ii]=posT[ii] - 2*pcell[1];
          }
          else if (posT[ii]<-pcell[1]){
              posT[ii]=posT[ii] + 2*pcell[1];
          }
      }

      for(ii=2;ii<3*Na;ii+=3){
          if (posT[ii]>pcell[2]){
              posT[ii]=posT[ii] - 2*pcell[2];
          }
          else if (posT[ii]<-pcell[2]){
              posT[ii]=posT[ii] + 2*pcell[2];
          }
      }

}
MKL_INT neighbours(int listchk, MKL_INT Na, double *posT, double *rcut, \
        double *pcell, double **rlistn, double **rn, MKL_INT rsize)
{
    double rcut2;
    double *rlist=*rlistn;
    double *r=*rn;
    rcut2=*rcut;
    rcut2=1.2*rcut2*rcut2;
//     printf("r1 = %.9f\n",r[0]);
//     printf("neigh %E\t%E\t%E\t%E\n",posT[0],posT[1],posT[2],posT[5]);
    if (listchk==1){
        MKL_INT i, ind;
        MKL_INT inc=1;
        MKL_INT *rcols;

          printf("before free neigh 1\n");
          mkl_free(rlist);
          printf("first free\n");
          mkl_free(r);

          printf("second free\n");

          rcols    = (MKL_INT *)mkl_calloc(Na, sizeof( MKL_INT ), 64);

          rlist    = (double *)mkl_calloc(1000*Na, sizeof( double ), 64);

          r    = (double *)mkl_calloc(1000*Na, sizeof( double ), 64);
  //         
          *rlistn=rlist;
          *rn=r;

          neigh1(Na,posT,rcut2,pcell,r,rlist,rcols);
  //         printf("r2 = %.9f\n",r[0]);
  //         printf(" ptr = %d\n",rsize);
  //         rsize=rcols[0];
  //         printf("rcol0 = %d\n",rcols[0]);
          ind=cblas_idamax (Na, rcols, inc);
  //         printf("ind = %d\n",ind);
          rsize=rcols[ind];
  //         printf("rsize = %d\n",rsize);
  //         for (i=0; i< Na; i++){
  //             if (rcols[i]>rsize){
  //                 rsize=rcols[i];
  //             }
  //         }

//         rlist=mkl_realloc(rlist,Na*rsize);
//         r=mkl_realloc(r,Na*rsize);

          mkl_free(rcols);
  //         mkl_free(rlist);
  //         mkl_free(r);
      }
      else if (listchk==0){
  //         r    = (double *)mkl_calloc((*rsize)*Na, sizeof( double ), 64);
          neigh2(Na,posT,rcut2,pcell,r,rlist,rsize);

//         mkl_free(r);
    }

      return rsize;

}
void dafem(int i,double *erate,double *pos,double *u, \
        double *v,double *a,double *pcell, int natoms,\
        double *rlist,double *r,int listchk,double *atomtype,\
        double *atomlocal,double *epsm,double *sigm,double *rcut, \
        int iterlim,double *pos0,double *dt,double *gamma, double *beta, \
        double *Vm,MKL_INT *Im,MKL_INT *Jm,double *un,double *vn, \
        double *an,double *Fc,double *Ec,double *Epot,double *rlistn, \
        double *rn,MKL_INT *rlistcolsn, double *mass, MKL_INT *iter, \
        double *dc, double *tole, double *tolf, double *told, int ntypes, double *posn, double *pcelln)
{
#define INFO   0
    MKL_INT     info=INFO, lda=natoms;
    MKL_INT     request=0;
    MKL_INT     sort=3;
    MKL_INT ii, Na, dof, dof2, locat, ibase1, ibase2;
    MKL_INT rsize;
    MKL_INT mtype=2; /*pardiso solver type for symmetric positive definite */
    double *posT;
    double *posT0;
    MKL_INT inc=1;
    MKL_INT job[8];
    MKL_INT NC[2];
    MKL_INT *AJm;
    MKL_INT *AJ;
    MKL_INT *AIm;
    MKL_INT *AI;
    double *Am;
    double *A;
    MKL_INT *AJk;
    MKL_INT *AIk;
    double *Ak;
    double *ddi;
    double *dd1;
    double *f1;
    double *residc;
    double *dis;
//     double *v2;
    void *pt[64];
    MKL_INT iparm[64]= {{0}};
    MKL_INT     maxfct=1;
    MKL_INT     mnum=1;
    MKL_INT perm;
    MKL_INT phase;
    MKL_INT nrhs = 1;
    MKL_INT msglvl = 1;
    MKL_INT error=0;
    double Ecn, Ecd, Fcn, Fcd1, Fcd2, Fcd, dcn, dcd, v2n, Ke, T;
    double kb = 8.6173324e-5;
    double *resid;

      Na=natoms;
      dof=3*Na;
      dof2=dof*dof;
  //     rsize=&rsize_space;
      rsize=0;
  //     printf("rsiz = %d\n",rsize);
  //     *rsize = 0;
      posT    = (double *)mkl_calloc(3*Na, sizeof( double ), 64);
      posT0    = (double *)mkl_calloc(3*Na, sizeof( double ), 64);
      f1    = (double *)mkl_calloc(3*Na, sizeof( double ), 64);
      residc    = (double *)mkl_calloc(3*Na, sizeof( double ), 64);
      dis    = (double *)mkl_calloc(3*Na, sizeof( double ), 64);
  //     v2    = (double *)mkl_calloc(3*Na, sizeof( double ), 64);
      resid    = (double *)mkl_calloc(3*Na, sizeof( double ), 64);

      cblas_dcopy(3*Na, pos, inc, posT, inc);
      cblas_dcopy(3*Na, pos0, inc, posT0, inc);
      char ordering, trans;
      double alpha = 1;

      ordering = 'c'; /*column major */
      trans = 't'; /* transpose*/

      mkl_dimatcopy (ordering, trans, Na, 3, alpha, posT, Na, 3);
      mkl_dimatcopy (ordering, trans, Na, 3, alpha, posT0, Na, 3);

//     i=i+1;
    double *ua; /*displacement of applied strain */

      ua    = (double *)mkl_calloc(3*natoms, sizeof( double ), 64);

      /*apply strain in the x direction */
      for(ii=0;ii<3*Na;ii+=3){
          ua[ii]=erate[0]*posT[ii];
      }
      /*apply strain in the y direction */
      for(ii=1;ii<3*Na;ii+=3){
          ua[ii]=erate[1]*posT[ii];
      }
      /*apply strain in the z direction */
      for(ii=2;ii<3*Na;ii+=3){
          ua[ii]=erate[2]*posT[ii];
      }

      vdAdd(3*Na,ua,u,u); /* u=ua+u; */

      vdAdd(3*Na,ua,posT,posT); /*D=D+ua; */

      /* Periodic cell straining */
      pcell[0]=pcell[0]+erate[0]*pcell[0];
      pcell[1]=pcell[1]+erate[1]*pcell[1];
      pcell[2]=pcell[2]+erate[2]*pcell[2];

      pbc(posT,pcell,Na);

      rsize=neighbours(listchk, Na, posT, rcut, pcell, &rlist, &r, rsize); 

      mkl_free(ua);
      mkl_free(posT);
      mkl_free(posT0);
      mkl_free(resid);
      mkl_free(u0);
      mkl_free(v0);
      mkl_free(a0);
      mkl_free(a1);
      mkl_free(di);
      mkl_free(dtemp);
      mkl_free(resid1);
  //     mkl_free(I);
  //     mkl_free(J);
  //     mkl_free(V);
  //     mkl_free(Ik);
  //     mkl_free(Jk);
  //     mkl_free(U_der2);
      mkl_free(Am);
      mkl_free(AJm);
      mkl_free(AIm);
  //     mkl_free(Ak);
  //     mkl_free(AJk);
  //     mkl_free(AIk);
  //     mkl_free(A);
  //     mkl_free(AJ);
  //     mkl_free(AI);
      mkl_free(ddi);
      mkl_free(dd1);
      mkl_free(f1);
      mkl_free(residc);
      mkl_free(dis);
      mkl_free(r);
      mkl_free(rlist);
  }
  void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
  {
      int natoms, i, listchk, iterlim, ntypes;
  //     size_t rlistcols;
      double *pos;
      double *pos0;
      double *dt;
      double *pcell;
      double *rlist;
      double *r;
      double *atomtype;
      double *atomlocal;
      double *epsm;
      double *sigm;
      double *rcut;
      double *Epot;
      MKL_INT *iter;
      MKL_INT *Im;
      MKL_INT *Jm;
      MKL_INT rlistcolsn;
      MKL_INT listsize;
      double *Vm;
      double *erate;
      double *u;
      double *v;
      double *a;
      double *mass;
      double *un;
      double *vn;
      double *an;
      double *gamma;
      double *beta;
      double *rn;
      double *rlistn;
      double *Ec;
      double *Fc;
      double *dc;
      double *tole;
      double *tolf;
      double *told;

      double *posn;
      double *pcelln;

     i=mxGetScalar(prhs[0]);
     erate=mxGetPr(prhs[1]);
     pos=mxGetPr(prhs[2]);
     u=mxGetPr(prhs[3]);
     v=mxGetPr(prhs[4]);
     a=mxGetPr(prhs[5]);
     pcell=mxGetPr(prhs[6]);
     natoms=mxGetScalar(prhs[7]);

     rlist=mxGetPr(prhs[8]);
     r=mxGetPr(prhs[9]);
     listchk=mxGetScalar(prhs[10]);
     atomtype=mxGetPr(prhs[11]);
     atomlocal=mxGetPr(prhs[12]);
     epsm=mxGetPr(prhs[13]);
     sigm=mxGetPr(prhs[14]);
     rcut=mxGetPr(prhs[15]);
     iterlim=mxGetScalar(prhs[16]);
     pos0=mxGetPr(prhs[17]);
     dt=mxGetPr(prhs[18]);
     gamma=mxGetPr(prhs[19]);
     beta=mxGetPr(prhs[20]);
     Vm=mxGetPr(prhs[21]);
     Im= (MKL_INT *)mxGetData(prhs[22]);
     Jm= (MKL_INT *)mxGetData(prhs[23]);
     mass=mxGetPr(prhs[24]);
     tole=mxGetPr(prhs[25]);
     tolf=mxGetPr(prhs[26]);

            

            
            
            

    

        

	
    

        
Best Answer
    

    

								

								
These lines show that rlist and r come from the MATLAB Memory Manager:
     rlist=mxGetPr(prhs[8]);
     r=mxGetPr(prhs[9]);
These values get passed to dafem, which then does this with them:
      mkl_free(r);
      mkl_free(rlist);
I am assuming the mkl Memory Manager is completely different from the MATLAB Memory Manager. So why would you expect this to work? And if mkl_free actually used the MATLAB Memory Manager, this would invalidate the data areas of some input variables which would lead to a crash.
Seems to me you need to, at the very least, comment these mkl_free lines out.
						    
						
		
    
        

    

Related Solutions
MATLAB: How to fill an mxArray in a MEX-file for MATLAB

								

								
This change has been incorporated into the documentation in Release 2010a (R2010a). For previous releases, read below for any additional information:
Below is code that shows four valid methods for filling an mxArray with data:
 
#include"mex.h"
#define ONE plhs[0]
#define TWO plhs[1]
#define THREE plhs[2]
#define FOUR plhs[3]
void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[] )
{
    const double data[]={2.1, 3.4, 2.3, 2.45};
    double *dynamicData1, *dynamicData2, *dynamicData4;
    double nonDynamicData[4];
    const int dims[]={0,0};
    int i;
/**********************************************************************************    
   METHOD 1, mxSetData with dynamic data
 *********************************************************************************/
      dynamicData1 = (double*)mxMalloc(4 * sizeof( double ));
      // Create the matrix as 0 by 0, since you are about to allocate the memory manually
      ONE = mxCreateNumericMatrix(0, 0, mxDOUBLE_CLASS, mxREAL);
      for( i = 0; i < 4; i++ )
      {
          // Here, you could be loading data from somewhere.
          // Instead, just copy from data[i]
          dynamicData1[i]=data[i];  
      }
      // Put the C array into the mxArray. You can do this because
      // the mxArray is empty (0x0) at this time, so there will be no
      // memory leak.
      mxSetData(ONE, dynamicData1);
      mxSetM( ONE, 2 );  
      mxSetN( ONE, 2 );
      // It is not necessary to call mxFree( dynamicData1 ) because 
      // that memory is pointed to inside the mxArray you are returning.
/**********************************************************************************    
   METHOD 2, mxSetPr with dynamic data, mxSetPr is only for double data 
 *********************************************************************************/
      dynamicData2 = (double*)mxMalloc(4 * sizeof( double ));
      // Create the matrix as 0 by 0, since you are about to allocate the 
      // memory manually.
      TWO = mxCreateNumericMatrix(0, 0, mxDOUBLE_CLASS, mxREAL);
      for( i = 0; i < 4; i++ )
      {
          // Here, you could be loading data from somewhere.
          // Instead, just copy from data[i]
          dynamicData2[i]=data[i];  
      }
      // Put the C array into the mxArray. You can do this because
      // the mxArray is empty (0x0) at this time so there will be no
      // memory leak
      mxSetPr( TWO, dynamicData2 );
      mxSetM( TWO, 2 );  
      mxSetN( TWO, 2 );
      // It is not necessary to call mxFree( dynamicData2 ) because 
      // that memory is pointed to inside the mxArray you are returning.
/**********************************************************************************    
   METHOD 3, mxGetPr with non-dynamic data
 *********************************************************************************/
      // Create the matrix as 2-by-2, since you are just copying existing data into it
      THREE = mxCreateNumericMatrix(2, 2, mxDOUBLE_CLASS, mxREAL);
      for( i = 0; i < 4; i++ )
      {
          // Here, you could be loading data from somewhere.
          // Instead, just copy from data[i]
          ((double*)mxGetPr(THREE))[i]=data[i];  
      }
/**********************************************************************************    
   METHOD 4, mxGetPr with dynamic data, using a copy of the data
 *********************************************************************************/
      dynamicData4 = (double*)mxMalloc(4 * sizeof( double ));
      // Create the matrix as 2-by-2, since you are just copying existing data into it
      FOUR = mxCreateNumericMatrix(2, 2, mxDOUBLE_CLASS, mxREAL);
      for( i = 0; i < 4; i++ )
      {
          // Here, you could be loading data from somewhere.
          // Instead, just copy from data[i]
          dynamicData4[i]=data[i];  
      }
      for( i = 0; i < 4; i++ )
      {
          // Copy the data in.
          ((double*)mxGetPr(FOUR))[i]=dynamicData4[i];  
      }
      // You must call mxFree( dynamicData4 ) at this time because 
      // that memory is not pointed to inside the mxArray you are 
      // returning. This is because you copied it into the array with mxGetPr.
      mxFree( dynamicData4 );
      return;
  }
The above example only uses the double data type, but you can also create mxArrays of other data type. The documentation has an example using unsigned 16-bit integers:
http://www.mathworks.com/help/releases/R2010a/techdoc/matlab_external/f12977.html#f19620
						    
						
MATLAB: How to output a 2D integer array using Matlab mex

								

								
I got the answer with the help from someone. I put it here for others' reference.
    // Output the results 
    plhs[0]                     = mxCreateNumericMatrix(col_outputBuff, row_outputBuff, mxINT32_CLASS, mxREAL);
    int* outputMatrix         = (int *)mxGetData(plhs[0]);
    // Read in the data
    for (int col=0; col < col_outputBuff; col++) {
        for (int row=0; row < row_outputBuff; row++) {
            outputMatrix[row + col*row_outputBuff] = outputBuff[col][row];
        }
    }
						    
						
        
        

            
            
            
Related Question
Can MEX BLAS library be used for native double matrix in C
Does the advanced C MEX-file example on using LAPACK and BLAS functions fail to compile and build correctly
Difficulties making wrapper mex file on linux
        
    

Answer 1

These lines show that rlist and r come from the MATLAB Memory Manager:

     rlist=mxGetPr(prhs[8]);
     r=mxGetPr(prhs[9]);

These values get passed to dafem, which then does this with them:

      mkl_free(r);
      mkl_free(rlist);

I am assuming the mkl Memory Manager is completely different from the MATLAB Memory Manager. So why would you expect this to work? And if mkl_free actually used the MATLAB Memory Manager, this would invalidate the data areas of some input variables which would lead to a crash.

Seems to me you need to, at the very least, comment these mkl_free lines out.