Matlab Mex C implementation of Cholesky Decomposition

Question

I am currently investigating runtimes of different matrix-inversion methods and therefore came across the Cholesky decomposition. In order to benchmark with the built-in cholesky-decomposition of matlab, I would like to transform my matlab-based implementation of the Cholesky decomposition into a C-Implementation with a Mex-Matlab-interface.

I tried my best with my limited programming skills in C and a bunch of tutorials and examples, but I can't get my Mex-interface to compile. I would really appreciate if someone could give me a hand.

Here is my code:

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

The Matlab-based implementation of Cholesky I am trying to implement in C is:

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

Many thanks!

Edit: Here is the fixed code if anyone is looking for a Matlab-mex-C-based implementation of the Cholesky Decomposition:

#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

First of all, it is not possible to multi-index an mxArray directly. In C, two dimensional matrices are arrays of arrays. You could create a macro or linearly compute the indices as both in this answer.

Index the matrices linearly: L[(i)+(k)*nrows)] instead of L[i][k].

The compiler tells you directly that you should

#include <math.h>

You also have to get the output pointer by

outMatrix = mxGetPr(plhs[0]);