Code:
double getlikelHood(double *r,double *sigtmp)
{
double *n;
int i,j;
double lwkopt;
int ctr;
int dim = DIMENSION;
double *work;
double alpha = 1;
double beta = 0;
double *y,*sigma;
char trans = 'N';
double y1;
double result;
double detVal,retVal;
int *iipiv,iinfo,wwork,iinc,mmone,lwork,*ipiv,size;
y = (double*)malloc(DIMENSION*sizeof(double));
sigma = (double*)malloc(DIMENSION*DIMENSION*sizeof(double));
n = (double*)malloc(DIMENSION*DIMENSION*sizeof(double));
mmone = -1;
iinc = 1;
for(i = 0;i<DIMENSION;i++)
{
for(j = 0;j<DIMENSION;j++)
{
*(sigma+i*DIMENSION+j) = *(sigtmp+i*DIMENSION+j);
*(n+i*DIMENSION+j) = *(sigtmp+i*DIMENSION+j);
}
}
detVal=0;
detVal=deter(n,DIMENSION);
ipiv = (int *)malloc(dim*dim*sizeof(int));
size = DIMENSION;
iinfo = 1;
dgetri_(&dim,n,&dim,ipiv,&lwkopt,&mmone,&iinfo);
lwork = (int)lwkopt;
work = (double*)malloc(lwork*lwork*sizeof(double));
dgetrf_(&dim,&dim,sigma,&size,ipiv,&iinfo);
dgetri_(&dim,sigma,&dim,ipiv,work,&lwork,&iinfo);
dgemv_(&trans,&size,&size,&alpha,sigma,&size,r,&iinc,&beta,y,&iinc);
result = -0.5*ddot_(&size,y,&iinc,r,&iinc);
retVal = exp(result)/sqrt(4*pi*pi*fabs(detVal));
free(y);
free(n);
free(work);
free(sigma);
free(ipiv);
return retVal;
}
].
