MKL library behaving differently in mex-files and in standalone C++

Question

I am trying to get GPc (https://github.com/SheffieldML/GPc) working in Matlab, using mex-files. I got the examples working, I took the bit I'm currently being interested in out as a standalone C++ program, that works just fine. However, when I try to do the same in a mex and run it through Matlab, I'm getting some errors, in particular:

MKL ERROR: Parameter 4 was incorrect on entry to DPOTRF.

or

** On entry to DPOTRF parameter number  4 had an illegal value

depending on whether I use the system version of MKL or the one Matlab carries along. The call to dpotrf is:

dpotrf_(type, nrows, vals, nrows, info);

with all variables valid (type="U", nrows=40, vals = double[40*40]) and with the interface:

extern "C" void dpotrf_(
        const char* t,  // whether upper or lower triangluar 'U' or 'L'
        const int &n,   // (input)
        double *a,  // a[n][lda] (input/output)
        const int &lda, // (input)
        int &info   // (output)
        );

(both are taken from GPc). LDA was originally supplied as ncols (which I believe is incorrect, but I didn't inquiry the library author about it yet), but it shouldn't make a difference, because this is called on a square matrix.

I feared that there might be problem with the references, so I changed the interface header to accept int* (like in http://www.netlib.org/clapack/clapack-3.2.1-CMAKE/SRC/dpotrf.c), but that started giving me segfaults, so it made me thinking the references there are right.

Does anybody have an idea what might be wrong?

score 2 · Accepted Answer · answered Dec 18 '14 at 11:39

I've tried to reproduce with an example on my end, but I'm not seeing any errors. In fact the result is identical to MATLAB's.

mex_chol.cpp

#include "mex.h"
#include "lapack.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    // verify arguments
    if (nrhs != 1 || nlhs > 1) {
        mexErrMsgTxt("Wrong number of arguments.");
    }
    if (!mxIsDouble(prhs[0]) || mxIsComplex(prhs[0])) {
        mexErrMsgTxt("Input must be a real double matrix.");
    }
    if (mxGetM(prhs[0]) != mxGetN(prhs[0])) {
        mexErrMsgTxt("Input must be a symmetric positive-definite matrix.");
    }

    // copy input matrix to output (its contents will be overwritten)
    plhs[0] = mxDuplicateArray(prhs[0]);

    // pointer to data
    double *A = mxGetPr(plhs[0]);
    mwSignedIndex n = mxGetN(plhs[0]);

    // perform matrix factorization
    mwSignedIndex info = 0;
    dpotrf("U", &n, A, &n, &info);

    // check if call was successful
    if (info < 0) {
        mexErrMsgTxt("Parameters had an illegal value.");
    } else if (info > 0) {
        mexErrMsgTxt("Matrix is not positive-definite.");
    }
}

Note that MATLAB already ships with BLAS/LAPCK headers and libraries (Intel MKL implementation). In fact this is what $MATLABROOT\extern\include\lapack.h has as function prototype for dpotrf:

#define dpotrf FORTRAN_WRAPPER(dpotrf)
extern void dpotrf(
    char   *uplo,
    ptrdiff_t *n,
    double *a,
    ptrdiff_t *lda,
    ptrdiff_t *info
);

Here is how you compile the above C++ code:

>> mex -largeArrayDims mex_chol.cpp libmwblas.lib libmwlapack.lib

Finally let's test the MEX function:

% some random symmetric positive semidefinite matrix
A = gallery('randcorr',10);

% our MEX-version of Cholesky decomposition
chol2 = @(A) triu(mex_chol(A));

% compare
norm(chol(A) - chol2(A))   % I get 0

(Note that the MEX code returns the working matrix as is, where the LAPACK routine only overwrites half of the matrix. So I used TRIU to zero-out the other half and extract the upper part).

I think you are confusing the Fortran interface with the C interface in your code — Amro, Dec 18 '14 at 11:49
@Amro Your comment is correct, I was compiling with a library using the C-interface (mkl_rt), but running with a library using Fortran interface (blas&lapack)! Thanks a lot! I would like to close the question somehow, but the correct answer is in the comment, so I can't really accept it like this, or how can I do it? — BIOStheZerg, Jan 09 '15 at 10:46
@Amro - Sorry to spam an old post, but I was curious to know if you know the answer to this question: http://stackoverflow.com/questions/29290232/quickly-and-efficiently-calculating-an-eigenvector-for-known-eigenvalue - It sounds like something up your alley and something you would know. — rayryeng, Mar 26 '15 at 23:29
@chappjc - Sorry to spam an old post, but I was curious to know if you know the answer to this question: http://stackoverflow.com/questions/29290232/quickly-and-efficiently-calculating-an-eigenvector-for-known-eigenvalue - It sounds like something up your alley and something you would know. — rayryeng, Mar 26 '15 at 23:30
@rayryeng Thanks for pointing out an interesting post. Luis to the rescue though! — chappjc, Mar 31 '15 at 04:33

MKL library behaving differently in mex-files and in standalone C++

1 Answers1

mex_chol.cpp