Prev: how to localize the pixel's position.
Next: PLP algorithm, few questions
From: abbas chatraei on 30 Apr 2010 06:55
Hi everybody, Please Help.

I wrote following C mex program. when I compile it, MATLAB dose not give any error. But when I call it from MATLAB prompt, MATLAB hangs.

I used some command from lapack library in this mex file and for compiling I used following commands:

lapacklib = fullfile(matlabroot, 'extern', 'lib', 'win32', 'microsoft', 'libmwlapack.lib');
mex ('-v', 'mylapck2.c', lapacklib);


call command:

at=0.2*ones(9,1);
val=[1 2 3 4 1];
[a,b,c,d]=mylapck2(at,val)



my c mex program:


*=======================================================*/
#include "mex.h"
#include "D:\Program Files\MATLAB\R2009a\extern\examples\refbook\fort.h"

void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )
{

/* mex interface to LAPACK functions dsysvx and zhesvx */

int m, N, n, lda, ldb, info;
double *at,*b, *a4, *U, *w, *x, *y, *z, *M, *MT, *kT, *bOut, *cOut, *dOut ;
double *diag1,*diag2,*diag3, *ones, *zeros, *zero1, *zero2, *k;
int i ,j, l=0;


at = mxGetPr(prhs[0]);
b = mxGetPr(prhs[1]);
m = mxGetM(prhs[0]);
n = mxGetN(prhs[0]);
N=n+1;



if (nrhs != 2) {
mexErrMsgTxt("Expect 2 input arguments and return 1 output argument");
}

U =(double *) mxCalloc((N-1)*N,sizeof(double));
w =(double *) mxCalloc(3*N+1,sizeof(double));
x =(double *) mxCalloc(3*N+1,sizeof(double));
y =(double *) mxCalloc(3*N+1,sizeof(double));
z =(double *) mxCalloc(3*N+1,sizeof(double));
M =(double *) mxCalloc((3*N+1)*(3*N+1),sizeof(double));
diag1 =(double *) mxCalloc((N-1)*(N-1), sizeof(double));
diag2 =(double *) mxCalloc((N-1)*(N-1), sizeof(double));
diag3 =(double *) mxCalloc((N-1)*(N-1), sizeof(double));
ones =(double *) mxCalloc(N-1, sizeof(double));
zeros =(double *) mxCalloc(N-1, sizeof(double));
zero1 =(double *) mxCalloc((N-1)*(N+1), sizeof(double));
zero2 =(double *) mxCalloc((N-1)*(2*N), sizeof(double));
k =(double *) mxCalloc((3*N+1), sizeof(double));





/********************************U matrix*/
for (i=0; i<N-1; i++){
for (j=0; j<N; j++){
U[i*(N)+j]=0;
}
}
for(i=0; i<N-1; i++){
U[i*(N)+i+1]=1;
U[i*(N)+i]=-1;
}
/**************************************/
for(i=0;i<(3*N+1); i++){
w[i]=0;
x[i]=0;
y[i]=0;
z[i]=0;
}
w[N-1]=b[4]*b[4]; w[2*N-1]=b[4]; w[3*N-1]=1; w[3*N]=b[4]*b[4]*b[4];
x[N-1]=2*b[4]; x[2*N-1]=1; x[3*N]=3*b[4]*b[4];
y[2*N]=1;
z[N]=1;
/*****************************************/
for(i=0;i<(N-1); i++){
ones[i]=0;
zeros[i]=0;
}
/*****************************************/
for (i=0; i<N-1; i++){
for (j=0; j<N-1; j++){
diag1[i*(N-1)+j]=0;
diag2[i*(N-1)+j]=0;
diag3[i*(N-1)+j]=0;
}
}
for(i=0;i<N-1;i++){
diag1[i*(N-1)+i]=-b[4]*b[4];
diag2[i*(N-1)+i]=-b[4];
diag3[i*(N-1)+i]=-2*b[4];
}
/********************************************/
for (i=0; i<N-1; i++){
for (j=0; j<N+1; j++){
zero1[i*(N+1)+j]=0;

}
}
for (i=0; i<N-1; i++){
for (j=0; j<2*N; j++){
zero2[i*(2*N)+j]=0;

}
}




/**[]N-1*N-1*****************/
for(i=0;i<N-1;i++)
for(j=0;j<N-1; j++){
M[i*(3*N)+j]=diag1[i*(N-1)+j];
M[i*(3*N)+N-1]=zeros[j];
M[i*(3*N)+j+N]=diag2[i*(N-1)+j];
M[i*(3*N)+2*N-1]=zeros[j];
M[i*(3*N)+j+2*N]=U[i*(N)+j];
M[i*(3*N)+3*N-1]=U[i*(N)+N-1];
M[i*(3*N)+3*N]=zeros[j];

M[(i+N-1)*(3*N)+j]=diag3[i*(N)+j];
M[(i+N-1)*(3*N)+N-1]=zeros[j];
M[(i+N-1)*(3*N)+j+N]=U[i*(N)+j];
M[(i+N-1)*(3*N)+2*N-1]=U[i*(N)+N-1];
M[(i+N-1)*(3*N)+j+2*N]=zero1[i*(N+1)+j];
M[(i+N-1)*(3*N)+3*N-1]=zero1[i*(N+1)+N-1];
M[(i+N-1)*(3*N)+3*N]=zero1[i*(N+1)+N];

M[(i+2*N-2)*(3*N)+j]=U[i*(N)+j];
M[(i+2*N-2)*(3*N)+ N-1]=U[i*(N)+N-1];
M[(i+2*N-2)*(3*N)+ N]=zeros[j];
M[(i+2*N-2)*(3*N)+ j+N+1]=zero2[i*(2*N)+j];
M[(i+2*N-2)*(3*N)+ j+2*N]=zero2[i*(2*N)+j+N-1];
M[(i+2*N-2)*(3*N)+ 3*N-1]=zero2[i*(2*N)+2*N-2];
M[(i+2*N-2)*(3*N)+ 3*N]=zero2[i*(2*N)+2*N-1];

M[(3*N-3)*(3*N)+ j]=w[j];

M[(3*N-2)*(3*N)+ j]=x[j];

M[(3*N-1)*(3*N)+ j]=y[j];

M[(3*N)*(3*N)+ j]=z[j];

}
/*****************************************/
for(i=0;i<(3*N+1); i++){
k[i]=0;
}
for(i=0;i<N-1;i++)
k[i]=at[i]*b[4]*b[4]*b[4];
for(i=N;i<2*(N-1);i++){
k[i]=3*at[l]*b[4]*b[4];
l++;
}

l=0;
for(i=2*(N-1);i<3*(N-1);i++){
k[i]=3*at[l]*b[4];
l++;
}
k[3*(N-1)+1]=b[2]; k[3*(N-1)+2]=b[3]; k[3*(N-1)+3]=b[0]; k[3*(N-1)+4]=b[1];
/********************************************************************************/
lda = 3*N+1;
ldb = 3*N+1;
MT = (double *)mxCalloc((3*N+1)*(3*N+1),sizeof(double));
kT = (double *)mxCalloc(3*N+1,sizeof(double));

for(i=0;i<3*N+1;i++) {
kT[i] = k[i];
}
for(i=0;i<(3*N+1)*(3*N+1);i++) {
MT[i] = M[i];
}

dgesv(&lda, &lda, MT, &lda, &lda, kT, &ldb, &info);

for(i=0;i<3*N+1;i++) {
mexPrintf("%f\t",kT[i]);
}

plhs[0] = mxCreateDoubleMatrix(N, 1, mxREAL);
plhs[1] = mxCreateDoubleMatrix(N, 1, mxREAL);
plhs[2] = mxCreateDoubleMatrix(N, 1, mxREAL);
plhs[3] = mxCreateDoubleMatrix(1, 1, mxREAL);
bOut = mxGetPr(plhs[0]);
cOut = mxGetPr(plhs[1]);
dOut = mxGetPr(plhs[2]);
a4 = mxGetPr(plhs[3]);

for(i=0;i<N;i++) {
bOut[i] = kT[i];
cOut[i] = kT[i+N-1];
dOut[i] = kT[i+2*N-1];
}
a4[0]= kT[3*N];
}
From: Daniel Armyr on 7 May 2010 03:23
For starters, you should have your program check your input agrumetns. For example, I see that you assume that b hasz at elast 5 elements. You shouldn't do that unless you have checked that first.

Over all, this is a rather messy piece of code. You may want to structure it a bit to help you out. For example, you have no idea where it crashes. Why can't you tell us on which line it crashes? Figure that out and you will have solved 90% of this problem.
 | 
Pages: 1
Prev: how to localize the pixel's position.
Next: PLP algorithm, few questions