Matlab Mex C Cholesky分解的实现

Question

我目前正在研究不同的矩阵求逆方法的运行时间，因此遇到了Cholesky分解。为了基于matlab的内置胆甾分解进行基准测试，我想将基于matlab的Cholesky分解实现转换为带有Mex-Matlab接口的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的实现是：

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 

非常感谢！

编辑：下面是如果有人正在寻找一个基于Matlab的MEX-C-执行Cholesky分解的固定代码：

#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); 
} 

Answer 1

首先，这是不可能的多直接索引mxArray。在C中，二维矩阵是arrays of arrays。你可以在this answer中创建宏或线性计算索引。

线性索引矩阵：L[(i)+(k)*nrows)]而不是L[i][k]。

编译器会告诉你直接说你应该

#include <math.h> 

您还可以通过

outMatrix = mxGetPr(plhs[0]);