Cholesky 分解的 Matlab Mex C 实现

Posted 2023-02-16

【中文标题】Cholesky 分解的 Matlab Mex C 实现

【英文标题】：Matlab Mex C implementation of Cholesky Decomposition

【发布时间】：2015-11-12 16:37:10

【问题描述】：

我目前正在研究不同矩阵求逆方法的运行时，因此遇到了 Cholesky 分解。为了与 matlab 的内置 cholesky 分解进行基准测试，我想将基于 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 的 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

非常感谢！

编辑：如果有人正在寻找基于 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);

【问题讨论】：

首先 L[i][i] = sqrt(M[i][i] - sum; 这里少了一个括号。

谢谢！我更正了它并添加了我的 Matlab 实现，以便更好地理解我试图实现的内容。

如果您发布编译器错误将非常有帮助。

抱歉，添加了编译器错误

好吧，我认为问题真的归结为矩阵运算，因为几乎每个错误都与此有关。就是不知道怎么引入矩阵，交给函数，对矩阵做运算，交还给主函数，然后返回。有什么想法吗？

【参考方案1】：

首先，不能直接对mxArray 进行多重索引。在 C 中，二维矩阵是arrays of arrays。您可以在this answer 中创建一个宏或线性计算索引。

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

编译器直接告诉你应该这样做

#include <math.h>

您还必须通过

获取输出指针

outMatrix = mxGetPr(plhs[0]);