Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
7117 Discussions

mex command for mkl linking on linux

Ali_R_3
Beginner
‎11-02-2016 11:59 AM
1,433 Views

Hi

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

 

0 Kudos

Link Copied

5 Replies
Ali_R_3
Beginner
‎11-02-2016 12:42 PM
1,433 Views

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

0 Kudos
Copy link
Gennady_F_Intel
Moderator
‎11-02-2016 08:55 PM
1,433 Views

>> 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

 

0 Kudos
Copy link
Ali_R_3
Beginner
‎11-03-2016 06:29 AM
1,433 Views

Hi

 

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. 

0 Kudos
Copy link
Ali_R_3
Beginner
‎11-03-2016 09:31 AM
1,433 Views

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)

 

0 Kudos
Copy link
Ying_H_Intel
Employee
‎11-08-2016 12:41 AM
1,433 Views

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 

0 Kudos
Copy link
Reply

Community support is provided Monday to Friday. Other contact methods are available here.

Intel does not verify all solutions, including but not limited to any file transfers that may appear in this community. Accordingly, Intel disclaims all express and implied warranties, including without limitation, the implied warranties of merchantability, fitness for a particular purpose, and non-infringement, as well as any warranty arising from course of performance, course of dealing, or usage in trade.

For more complete information about compiler optimizations, see our Optimization Notice.