mexGetPr error when freeing the data

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]);
     told=mxGetPr(prhs[27]);
     ntypes=mxGetScalar(prhs[28]);

//     rlistcols = mxGetN(prhs[8]);

      plhs[0] = mxCreateDoubleMatrix(3*natoms,1, mxREAL);
      plhs[1] = mxCreateDoubleMatrix(3*natoms,1, mxREAL);
      plhs[2] = mxCreateDoubleMatrix(3*natoms,1, mxREAL);
      plhs[3] = mxCreateDoubleMatrix(1,1, mxREAL);
      plhs[4] = mxCreateDoubleMatrix(1,1, mxREAL);
      plhs[5] = mxCreateDoubleMatrix(1,1, mxREAL);
      plhs[9] = mxCreateDoubleMatrix(1,1, mxREAL);

      plhs[10] = mxCreateDoubleMatrix(natoms,3, mxREAL);
      plhs[11] = mxCreateDoubleMatrix(3,1, mxREAL);
      posn = mxGetPr(plhs[10]);
      pcelln = mxGetPr(plhs[11]);

//     plhs[10] = mxCreateDoubleMatrix(natoms,3, mxREAL);
//     plhs[11] = mxCreateDoubleMatrix(3,1, mxREAL);

      rlistn = mxCalloc(natoms*1000, sizeof(double)); /* may have to change
                                                        to MKL_INT*/
      plhs[6] = mxCreateDoubleMatrix(0,0, mxREAL);

      rn = mxCalloc(natoms*1000, sizeof(double));
      plhs[7] = mxCreateDoubleMatrix(0,0, mxREAL);

      iter = mxCalloc(1, sizeof(MKL_INT));
      plhs[8] = mxCreateNumericMatrix(0, 0, mxUINT64_CLASS, mxREAL);
      mxSetData(plhs[8], iter);
      mxSetM(plhs[8], 1);
      mxSetN(plhs[8], 1);

      un = mxGetPr(plhs[0]);
      vn = mxGetPr(plhs[1]);
      an = mxGetPr(plhs[2]);
      Fc = mxGetPr(plhs[3]);
      Ec = mxGetPr(plhs[4]);
      Epot = mxGetPr(plhs[5]);
      dc = mxGetPr(plhs[9]);

//     posn=mxGetPr(plhs[10]);
//     pcelln=mxGetPr(plhs[11]);

//     plhs[10] = mxCreateDoubleMatrix(natoms,3, mxREAL);
//     plhs[11] = mxCreateDoubleMatrix(3,1, mxREAL);

//     *rlistcolsn=1000;

      dafem(i,erate,pos,u,v,a,pcell,natoms,\
          rlist,r,listchk,atomtype,atomlocal,epsm,sigm,rcut,iterlim, \
          pos0,dt,gamma,beta,Vm,Im,Jm,un,vn,an,Fc,Ec, \
              Epot,rlistn,rn,&rlistcolsn, mass, iter, dc, tole, tolf, told, ntypes, posn, pcelln);

      listsize=natoms*rlistcolsn;
      rlistn = (double*)mxRealloc(rlistn, sizeof(double)*listsize);

      mxSetData(plhs[6], rlistn);
      mxSetM(plhs[6], natoms);
      mxSetN(plhs[6], rlistcolsn);

      rn = (double*)mxRealloc(rn, sizeof(double)*listsize);

      mxSetData(plhs[7], rn);
      mxSetM(plhs[7], natoms);
      mxSetN(plhs[7], rlistcolsn);
  }

Answer 1

[James Tursa]

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.

[Ali Radhi]

Hi

Thank you for the the quick reply. I am a beginner at this. I already know that this won't work. But my question is that why does it work on windows but not linux. Also, i need to free the data at first so i can't comment those lines out. The only way i managed to make my code work is to add the following before the free statement:

r=NULL; rlist=NULL;

Could you explain?

[James Tursa]

Memory managers on linux work differently than memory managers on Windows. The fact that one system responds to the incorrect code by crashing while the other system doesn't does not mean anything ... it is still incorrect code in both cases. You just got lucky on one system that it didn't crash right away (could be that something is still messed up in memory and it will crash later).

And doing this:

r=NULL; rlist=NULL;

prior to calling mkl_free is equivalent to commenting the mkl_free lines out. (i.e., calling mkl_free(NULL) does nothing)

[Ali Radhi]

Ok. I understand. One last question then. The reason I do the free is to create a new r and rlist variables with neigh1 call in neighbours function, while destroying the old one. That is why I do calloc afterwards immediately. If I comment out the mkl_free statements, will calloc delete the old one and create a new one at the same time? I only ask because I do not want to have increased memory (memory leak) every time I call this function. I thank you for your help. As you know, I am a bit of a beginner at this. Thank you again.

[James Tursa]

If r and rlist are sometimes coming from a prhs[] variable and sometimes coming from mkl_calloc, then to manage this correctly you will need to keep track of where it came from in your code. If it came from a prhs[] variable, then don't free it! But if it came from a mkl_calloc call, then free it. The calloc function will not free any existing allocated memory.

[Ali Radhi]

Hi. Thanks for the reply. If my variables came from prhs[], then how to free its space?

[James Tursa]

The prhs[ ] variables are intended to be READ-ONLY, hence the const qualifier. You are not supposed to free its memory inside a mex routine because you can screw up other variables that might be sharing data with this variable, leading to a crash.