topic I found out a segmentation in Intel® oneAPI Math Kernel Library

mex command for mkl linking on linux

Ali_R_3 — Wed, 02 Nov 2016 18:59:38 GMT

I am trying to compile a c code through mex function in matlab. I was able to do so on a windows machine, but not on a linux one. The mkl line advisor gave me the following:

Link line:

-Wl,--start-group ${MKLROOT}/lib/intel64/libmkl_intel_ilp64.a ${MKLROOT}/lib/intel64/libmkl_intel_thread.a ${MKLROOT}/lib/intel64/libmkl_core.a -Wl,--end-group -liomp5 -lpthread -lm -ldl

Compiler options:

-DMKL_ILP64 -m64 -I${MKLROOT}/include

I tried many things, but nothing helped. Followingsome references in the net, i used the following makefile:

#!/bin/make

BINARY = mexDAFEM.mexa64

MATLAB = /usr/local/packages/matlabR2016a
MKLROOT = /u/c/radhiali/intel/mkl
MKLLIB = $(MKLROOT)/lib/intel64

CC = gcc
OPTS = -DMATLAB_MEX_FILE -DMKL_ILP64 -fPIC -fno-omit-frame-pointer -DMX_COMPAT_32 -O3 -D_GNU_SOURCE
SRC = mexDAFEM.c
OBJS = $(SRC:.c=.o)

INC = -I$(MATLAB)/extern/include -I$(MATLAB)/simulink/include -I$(MKLROOT)/include
LIBS = -shared -Wl,--version-script,$(MATLAB)/extern/lib/glnxa64/mexFunction.map -Wl,--no-undefined
LIBS += -Wl,-rpath-link,$(MATLAB)/bin/glnxa64 -L$(MATLAB)/bin/glnxa64 -lmx -lmex -lmat
LIBS += -Wl,--start-group $(MKLLIB)/libmkl_intel_ilp64.a $(MKLLIB)/libmkl_sequential.a $(MKLLIB)/libmkl_core.a -Wl,--end-group
LIBS += -lpthread -lm -ldl

.SUFFIXES: .c .o

.c.o:
   $(CC) $(OPTS) $(INC) -c $< -o $@

$(BINARY): $(OBJS)
   $(CC) $(OBJS) $(OPTS) -o $(BINARY) $(LIBS)
   @echo Binary created!!

It compiled, but the matlab shows a system error, saying it encountered an internal problem and needs to close.

Please advise as I am at my wit's end with this problem. Thank you in advance for making the time.

Regards

Ali Radhi

I also wanted to mention that

Ali_R_3 — Wed, 02 Nov 2016 19:42:35 GMT

I also wanted to mention that My compiler is

gcc (GCC) 4.4.7 20120313 (Red Hat 4.4.7-17)

Matlab is 2016a.

The code I wanted to compile is as follows:

#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;
        
        
        mkl_free(rlist);
        mkl_free(r);
        
      
        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>rsize){
//                 rsize=rcols;
//             }
//         }
        
//         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 main(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); 
    
    /* force calculation */

    
    *Epot=nbresidual(Na,ntypes,0,posT,rlist,r,atomtype,atomlocal,pcell,epsm,sigm,rcut,resid,rsize);
    
    double *u0;
    double *v0;
    double *a0;
    double *a1;
    double *di;
    double *dtemp;
    double *resid1;

    u0    = (double *)mkl_calloc(3*Na, sizeof( double ), 64);
    v0    = (double *)mkl_calloc(3*Na, sizeof( double ), 64);
    a0    = (double *)mkl_calloc(3*Na, sizeof( double ), 64);
    di    = (double *)mkl_calloc(3*Na, sizeof( double ), 64);
    dtemp    = (double *)mkl_calloc(3*Na, sizeof( double ), 64);
    resid1    = (double *)mkl_calloc(3*Na, sizeof( double ), 64);

    a1    = (double *)mkl_calloc(3*Na, sizeof( double ), 64);
    
    cblas_dcopy(3*Na, resid, inc, a, inc);
    
    double mi;
    
    mi=1/(*mass);
    
 
    trans = 'n'; /* normal*/
    
//     mkl_dimatcopy (ordering, trans, 3*Na, 1, mi, a, 3*Na, 3*Na);
    cblas_dscal (dof, mi, a, inc);
    
    cblas_dcopy(3*Na, u, inc, u0, inc);
    cblas_dcopy(3*Na, v, inc, v0, inc);
    cblas_dcopy(3*Na, a, inc, a0, inc);
    cblas_dcopy(3*Na, u, inc, di, inc);
    
    *iter=0;
    /* start of minimization loop */
    MKL_INT chk=0; /* breaking loop flag */
    
    MKL_INT ntriplets;
    MKL_INT *Ik;
    MKL_INT *Jk;
    double *U_der2;
//     MKL_INT *I;
//     MKL_INT *J;
//     double *V;
    
    locat=2;/* for coo to csr*/
    ibase1=1;
    ibase2=1;
    job[1]=ibase1;
    job[2]=ibase2;
    job[3]=locat;
    job[4]=dof;
    job[0]=2;
    job[5]=0;
    
    AJm    = (MKL_INT *)mkl_calloc(dof, sizeof( MKL_INT ), 64);
    AIm    = (MKL_INT *)mkl_calloc(dof+1, sizeof( MKL_INT ), 64);
    Am    = (double *)mkl_calloc(dof, sizeof( double ), 64);
    
    mkl_dcsrcoo(job, &dof, Am, AJm, AIm, &dof, Vm, Im, Jm, &info);
    
    double betafac, a0fac;
    
//     betafac = *dt;
//     betafac = *beta*betafac*betafac;
//     betafac = 1/betafac;
    betafac=1/((*beta)*(*dt)*(*dt));
    
//     a0fac = *dt;
//     a0fac = a0fac*a0fac*(-0.5+*beta);
    a0fac = -0.5*(*dt)*(*dt)*(1-2*(*beta));
    
    double gammafac1, gammafac2;
//     gammafac1 = *dt;
//     gammafac1 = gammafac1*(1-*gamma);
    gammafac1 = (1-*gamma)*(*dt);
//     gammafac2 = *dt;
//     gammafac2 = gammafac2*(*gamma);
    gammafac2 = (*gamma)*(*dt);

    
    ddi    = (double *)mkl_calloc(dof, sizeof( double ), 64);
    dd1    = (double *)mkl_calloc(dof, sizeof( double ), 64);
    
    pardisoinit (pt,&mtype,iparm);
    iparm[5] = 0;
    
    while (*iter < iterlim && chk == 0)
    {
        
        *iter=*iter+1;
//         vdAdd(3*Na,di,posT0,posT); /* D=D0+di; */
//         pbc(posT,pcell,Na);
        
        rsize=neighbours(listchk, Na, posT, rcut, pcell, &rlist, &r, rsize);
        
        *Epot=nbresidual(Na,ntypes,0,posT,rlist,r,atomtype,atomlocal,pcell,\
                epsm,sigm,rcut,resid,rsize);
        

        
        /* a=(1/(beta*dt*dt))*(di-u0-dt*v0-0.5*dt*dt*(1-2*beta)*a0); */
        vdSub( 3*Na, di, u0, dtemp );
        cblas_daxpy (3*Na, -(*dt), v0, inc, dtemp, inc);
        cblas_daxpy (3*Na, a0fac, a0, inc, dtemp, inc);
//         mkl_dimatcopy (ordering, trans, 3*Na, 1, betafac, dtemp,3*Na,3*Na);
        cblas_dscal (dof, betafac, dtemp, inc);
        cblas_dcopy(3*Na, dtemp, inc, a, inc);
        
        /* v=v0+(1-gamma)*dt*a0+gamma*dt*a; */
        
        cblas_dcopy(3*Na, v0, inc, dtemp, inc);
        
        cblas_daxpy (3*Na, gammafac1, a0, inc, dtemp, inc);
        cblas_daxpy (3*Na, gammafac2, a, inc, dtemp, inc);
        
        cblas_dcopy(3*Na, dtemp, inc, v, inc);
        
        /* Ri=fa-M*a; */

        cblas_daxpy (3*Na, -(*mass), a, inc, resid, inc); /*resid = Ri */
        
        /* Stiffness matrix Calculation */   
        ntriplets=0;
        
        Ik    = (MKL_INT *)mkl_calloc(50e6, sizeof( MKL_INT ), 64);
        Jk    = (MKL_INT *)mkl_calloc(50e6, sizeof( MKL_INT ), 64);
        U_der2    = (double *)mkl_calloc(50e6, sizeof( double ), 64);        

//         
        ntriplets=nbatomstiffness(Na,rlist,r,posT,atomtype,atomlocal,pcell,epsm,\
                sigm,rcut,Ik, Jk, U_der2, rsize);
        
//         nbatomstiffness(Na,rlist,r,posT,atomtype,atomlocal,pcell,epsm,\
//                 sigm,rcut,I, J, V,ntriplets, rsize);

        
//         I=mkl_realloc(I,*ntriplets);
//         J=mkl_realloc(J,*ntriplets);
//         V=mkl_realloc(V,*ntriplets);
        
        /* from Coo to CSR sparse formating for stiffness*/
        job[4]=ntriplets;
        
        AJk    = (MKL_INT *)mkl_calloc(dof*dof, sizeof( MKL_INT ), 64);
        AIk    = (MKL_INT *)mkl_calloc(dof+1, sizeof( MKL_INT ), 64);
        Ak    = (double *)mkl_calloc(dof*dof, sizeof( double ), 64);
        
        mkl_dcsrcoo(job, &dof, Ak, AJk, AIk, &ntriplets, U_der2, Ik, Jk, &info);
        
        mkl_free(Ik);
        mkl_free(Jk);
        mkl_free(U_der2);
        
        NC[0]=0;
//         
//         for (ii=0; ii<ntriplets; ii++){
//             if (abs(AJk[ii])>1.0e-100){
//                 NC[0]=NC[0]+1;
//             }
//         }
//         
//         Ak=mkl_realloc(Ak,NC[0]);
//         AJk=mkl_realloc(AJk,NC[0]);
        
        /* A=(1/(beta*dt*dt))*M+Ka; */
        AJ    = (MKL_INT *)mkl_calloc(dof*dof, sizeof( MKL_INT ), 64);
        AI    = (MKL_INT *)mkl_calloc(dof+1, sizeof( MKL_INT ), 64);
        A    = (double *)mkl_calloc(dof*dof, sizeof( double ), 64);

        
        mkl_dcsradd (&trans,&request,&sort,&dof,&dof,Ak,AJk,AIk,&betafac\
                ,Am,AJm,AIm,A,AJ,AI,&ntriplets,&info);
        
        mkl_free(Ak);
        mkl_free(AJk);
        mkl_free(AIk);
        NC[1]=0;
        
//         for (ii=0; ii<ntriplets; ii++){
//             if (abs(AJ[ii])>1.0e-100){
//                 NC[1]=NC[1]+1;
//             }
//         }
//         A=mkl_realloc(Ak,NC[1]);
//         AJ=mkl_realloc(AJk,NC[1]);
        
//         A=mkl_realloc(A,NC[1]);
//         AJ=mkl_realloc(AJ,NC[1]);
//         
                phase=11;
//         printf("maxfct %d mnum %d mtype %d phase %d dof %d perm %d nrhs %d msglvl %d\n",maxfct\
//                 ,mnum,mtype,phase,dof,perm,nrhs,msglvl);
        pardiso_64 (pt,&maxfct,&mnum,&mtype,&phase,&dof,A,AI,AJ,&perm,&nrhs ,iparm,&msglvl,resid,ddi,&error);
        
//         printf("Perror1 = %d\n",error);
        
        phase=22;
        
        pardiso_64 (pt,&maxfct,&mnum,&mtype,&phase,&dof,A,AI,AJ,&perm,&nrhs ,iparm,&msglvl,resid,ddi,&error);
        
//         printf("Perror2 = %d\n",error);
        
        phase=33;
        
        pardiso_64 (pt,&maxfct,&mnum,&mtype,&phase,&dof,A,AI,AJ,&perm,&nrhs ,iparm,&msglvl,resid,ddi,&error);
        
//         printf("Perror3 = %d\n",error);
        

            phase=-1;
    
    pardiso_64 (pt,&maxfct,&mnum,&mtype,&phase,&dof,A,AI,AJ,&perm,&nrhs ,iparm,&msglvl,resid,ddi,&error);

        mkl_free(A);
        mkl_free(AJ);
        mkl_free(AI);

        
//         printf("Perror4 = %d\n",error);
        
        vdAdd(dof,ddi,di,di); /* di=ddi+di; */
        
        vdAdd(3*Na,di,posT0,posT); /* D=D0+di; */
        pbc(posT,pcell,Na);
        
        rsize=neighbours(listchk, Na, posT, rcut, pcell, &rlist, &r, rsize);
        
        *Epot=nbresidual(Na,ntypes,0,posT,rlist,r,atomtype,atomlocal,pcell,\
                epsm,sigm,rcut,residc,rsize);
        
        
        if (*iter == 1){
            
            cblas_dcopy(3*Na, resid, inc, resid1, inc); /* R1=Ri; */
            cblas_dcopy(3*Na, ddi, inc, dd1, inc); /* dd1=ddi; */
            cblas_dcopy(3*Na, a, inc, a1, inc); /* a1=M*a;*/
            cblas_dscal (dof, *mass, a1, inc);
//             cblas_daxpy (3*Na, *mass, a, inc, a1, inc); /* a1=M*a;*/
            cblas_dcopy(3*Na, residc, inc, f1, inc); /*f1=fc;*/
            
            Ecd = cblas_ddot (dof, dd1, inc, resid1, inc);
            printf(" ECD = %.9f\n",Ecd);
            Fcd1 = cblas_dnrm2 (dof, a1, inc);
            Fcd2 = cblas_dnrm2 (dof, f1, inc);
            if (Fcd1 > Fcd2){
                Fcd=Fcd1;
            }
            else {
                Fcd=Fcd2;
            }
            
        }
        
        
        cblas_daxpy (3*Na, -(*mass), a, inc, residc, inc); /* Rc=fc-M*a; */
        
        Ecn = cblas_ddot (dof, ddi, inc, residc, inc);
        
        *Ec = Ecn/Ecd;  /* Ec=(ddi'*Rc)/(dd1'*R1); */
        
        Fcn = cblas_dnrm2 (dof, residc, inc);
        
        *Fc = Fcn/Fcd; /* Fc=norm(Rc)/max(norm(M*a1),norm(f1)); */
        
        vdAdd(3*Na,di,u0,dis);
        
        dcn = cblas_dnrm2 (dof, ddi, inc);
        dcd = cblas_dnrm2 (dof, dis, inc);
        
        *dc = dcn/dcd;
        
        
//         vdSqr( dof, v, v2 );
        v2n = cblas_dnrm2 (dof, v, inc);
        v2n=v2n*v2n;
        
        Ke = 0.5*(*mass)*v2n; /* Ke = 0.5 v'*M*v only works of single species*/
        
        T= (2*Ke)/(kb*(dof-3));
        
        
        /*printf(out,'%d\t%d\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n',
         * i,iter,di(1),dc,Fc,Ec,Epot,0.5*v'*M*v,mass*sum(v.^2)/(kb*(3*natoms-3))); */
       printf("%d\t%d\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n" \
               ,i,*iter,di[0],*dc,*Fc,*Ec,*Epot,Ke,T);
       if (*Ec <=*tole && *Fc<=*tolf && *dc<=*told){
           
           cblas_dcopy(3*Na, di, inc, u, inc); /* u=di; */
           chk=1;                              /* break */
       }

    }
    
//     phase=-1;
    
//     pardiso_64 (pt,&maxfct,&mnum,&mtype,&phase,&dof,A,AI,AJ,&perm,&nrhs ,iparm,&msglvl,resid,ddi,&error);
    
    trans = 't'; /* transpose*/
    
    mkl_dimatcopy (ordering, trans, 3, Na, alpha, posT, 3, Na);
//     mkl_dimatcopy (ordering, trans, 3, Na, alpha, posT0, 3, Na);
    
//     cblas_dcopy(3*Na, posT, inc, pos, inc);
//     cblas_dcopy(3*Na, posT0, inc, pos0, inc);
    
    cblas_dcopy(3*Na, posT, inc, posn, inc);
    cblas_dcopy(3, pcell, inc, pcelln, inc);
    
    cblas_dcopy(3*Na, a, inc, an, inc);
    cblas_dcopy(3*Na, v, inc, vn, inc);
    cblas_dcopy(3*Na, u, inc, un, inc);
    
    cblas_dcopy(rsize*Na, r, inc, rn, inc);
    cblas_dcopy(rsize*Na, rlist, inc, rlistn, inc);
    *rlistcolsn=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;
    
    main(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);
}

Thank you.

Regards

Ali Radhi

>> It compiled, but the

Gennady_F_Intel — Thu, 03 Nov 2016 03:55:11 GMT

>> It compiled, but the matlab shows a system error, saying it encountered an internal problem and needs to close...

Ali, You never check error status when call MKL's functions. Could you add this and check which mkl's function produces the problem? or you may try to comment some of mkl's routines and check how it will work ......

You may also try to create the standalone C example and check if the problem is reproduced with the latest MKL ( vv2017 or 11.3.4 ). AFAIK R2016a uses MKL 11.1.1, but may be I am wrong. you may check by calling "version -blas" command from Matlab,

regards

Hi

Ali_R_3 — Thu, 03 Nov 2016 13:29:54 GMT

I am trying to debug step by step. I found out that the mkl_free function in neighbours functions in line 66 and 67 gives me an error. I tried commenting but the error goes on into the neigh1 function, which is a simple two embedded loops. When the first loop tries to go over the seconds, it gives an error. Any thoughts.

I found out a segmentation

Ali_R_3 — Thu, 03 Nov 2016 16:31:07 GMT

I found out a segmentation fault in line 65 or 66 (mk_free function) . any thoughts? My latest makefile looks like this:

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)

Hi Ali,

Ying_H_Intel — Tue, 08 Nov 2016 08:41:55 GMT

Hi Ali,

Do you want to Matlab and MKL 's integer type as 64bit? ( seems no -largeArrayDims is selected)

How about try the LP64 mode (without -DMKL_ILP64)

CPPFLAGS = -DMX_COMPAT_32  -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_lp64.a \
           $(MKL)/lib/intel64/libmkl_sequential.a \
           $(MKL)/lib/intel64/libmkl_core.a \
           -Wl,--end-group

Best Regards,
Ying