我目前正在研究不同矩阵求逆方法的运行时间,因此遇到了 Cholesky 分解。为了与 matlab 的内置 cholesky 分解进行基准测试,我想将基于 matlab 的 Cholesky 分解实现转换为具有 Mex-Matlab 接口(interface)的 C 实现。
我用我有限的 C 编程技能和一堆教程和示例尽了最大努力,但我无法编译我的 Mex 界面。如果有人能帮助我,我将不胜感激。
这是我的代码:
#include "mex.h"
/* The computational routine */
void myCholeskyC(double *M, double *L, mwSize n)
{
mwSize i;
mwSize j;
mwSize k;
mwSize l;
/* multiply each element y by x */
for (i=0; i<n; i++) {
double sum = 0;
for (k=0; k<n; k++) {
sum+= L[i][k]*L[i][k];
} //end of for-loop k
L[i][i] = sqrt(M[i][i] - sum);
for (j=i+1; j<n; j++) {
double sum2 = 0;
for (l=0; l<n; l++) {
sum2+= L[i][l]*L[j][l];
} //end of for-loop l
L[j][i] = (M[j][i] - sum2)/L[i][i];
} //end of for-loop j
} //end of for-loop i
} //end of computational routine
/* The gateway function */
void mexFunction( int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
{
double multiplier; /* input scalar */
double *inMatrix; /* 1xN input matrix */
size_t ncols; /* size of matrix */
double *outMatrix; /* output matrix */
/* check for proper number of arguments */
if(nrhs!=1) {
mexErrMsgIdAndTxt("MyToolbox:arrayProduct:nrhs","Only one input required.");
}
if(nlhs!=1) {
mexErrMsgIdAndTxt("MyToolbox:arrayProduct:nlhs","One output required.");
}
/* make sure the input argument is type double */
if( !mxIsDouble(prhs[0]) ||
mxIsComplex(prhs[0])) {
mexErrMsgIdAndTxt("MyToolbox:arrayProduct:notDouble","Input matrix must be type double.");
}
/* create a pointer to the real data in the input matrix */
inMatrix = mxGetPr(prhs[0]);
/* get dimensions of the input matrix */
ncols = mxGetN(prhs[0]);
/* create the output matrix */
plhs[0] = mxCreateDoubleMatrix((mwSize)ncols,(mwSize)ncols,mxREAL);
/* set output matrix */
outMatrix = plhs[0];
/* call the computational routine */
myCholeskyC(inMatrix,outMatrix,(mwSize)ncols);
}
我试图在 C 中实现的基于 Matlab 的 Cholesky 实现是:
function L = cholMatlab2(M)
n = length( M );
L = zeros( n, n );
for i=1:n
L(i, i) = sqrt(M(i, i) - L(i, :)*L(i, :)');
for j=(i + 1):n
L(j, i) = (M(j, i) - L(i,:)*L(j ,:)')/L(i, i);
end
end
end
非常感谢!
编辑:如果有人正在寻找 Cholesky 分解的基于 Matlab-mex-C 的实现,这里是固定代码:
#include "mex.h"
#include <math.h>
#define L(x,y) L[(x) + (y)*n]
#define M(x,y) M[(x) + (y)*n]
/* The computational routine */
void myCholeskyC(double *M, double *L, mwSize n)
{
mwSize i;
mwSize j;
mwSize k;
mwSize l;
/* multiply each element y by x */
for (i=0; i<n; i++) {
double sum = 0;
for (k=0; k<n; k++) {
sum+= L(i,k)*L(i,k);
} //end of for-loop k
L(i,i) = sqrt(M(i,i) - sum);
for (j=i+1; j<n; j++) {
double sum2 = 0;
for (l=0; l<n; l++) {
sum2+= L(i,l)*L(j,l);
} //end of for-loop l
L(j,i) = (M(j,i) - sum2)/L(i,i);
} //end of for-loop j
} //end of for-loop i
} //end of computational routine
/* The gateway function */
void mexFunction( int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
{
double multiplier; /* input scalar */
double *inMatrix; /* 1xN input matrix */
size_t ncols; /* size of matrix */
double *outMatrix; /* output matrix */
/* check for proper number of arguments */
if(nrhs!=1) {
mexErrMsgIdAndTxt("MyToolbox:arrayProduct:nrhs","Only one input required.");
}
if(nlhs!=1) {
mexErrMsgIdAndTxt("MyToolbox:arrayProduct:nlhs","One output required.");
}
/* make sure the input argument is type double */
if( !mxIsDouble(prhs[0]) ||
mxIsComplex(prhs[0])) {
mexErrMsgIdAndTxt("MyToolbox:arrayProduct:notDouble","Input matrix must be type double.");
}
/* create a pointer to the real data in the input matrix */
inMatrix = mxGetPr(prhs[0]);
/* get dimensions of the input matrix */
ncols = mxGetN(prhs[0]);
/* create the output matrix */
plhs[0] = mxCreateDoubleMatrix((mwSize)ncols,(mwSize)ncols,mxREAL);
/* set output matrix */
outMatrix = mxGetPr(plhs[0]);
/* call the computational routine */
myCholeskyC(inMatrix,outMatrix,(mwSize)ncols);
}
最佳答案
首先,无法直接对mxArray
进行多重索引。在 C 中,二维矩阵是 arrays of arrays .您可以在 this answer 中创建宏或线性计算索引.
线性索引矩阵:L[(i)+(k)*nrows)]
而不是 L[i][k]
。
编译器直接告诉你应该
#include <math.h>
您还必须通过以下方式获取输出指针
outMatrix = mxGetPr(plhs[0]);
关于c - Cholesky 分解的 Matlab Mex C 实现,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/33676622/