MKL 库在 mex 文件和独立C++中的行为不同

Question

我正在尝试使用mex文件让GPc(https://github.com/SheffieldML/GPc(在Matlab中工作。我让示例正常工作，我将我目前感兴趣的部分作为独立的C++程序，效果很好。但是，当我尝试在 mex 中执行相同的操作并通过 Matlab 运行它时，我遇到了一些错误，特别是：

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

或

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

取决于我使用的是 MKL 的系统版本还是 Matlab 附带的版本。对dpotrf的调用是：

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

所有变量都有效(type="U"，nrows=40，vals = double[40*40](，接口：

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

(两者都取自GPc(。LDA 最初是作为 ncols 提供的(我认为这是不正确的，但我还没有询问库作者(，但它应该没有区别，因为这是在方阵上调用的。

我担心引用可能有问题，所以我更改了接口标头以接受 int*(就像在 http://www.netlib.org/clapack/clapack-3.2.1-CMAKE/SRC/dpotrf.c 中一样(，但这开始给我带来段错误，所以它让我认为那里的引用是正确的。

有人知道可能出了什么问题吗？

Answer 1

我试图用一个例子来重现，但我没有看到任何错误。事实上，结果与MATLAB的结果相同。

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.");
    }
}

请注意，MATLAB 已经附带了 BLAS/LAPCK 标头和库(英特尔 MKL 实施(。事实上，这就是$MATLABROOTexternincludelapack.h作为dpotrf的功能原型：

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

以下是编译上述C++代码的方法：

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

最后，让我们测试 MEX 函数：

% 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

(请注意，MEX 代码按原样返回工作矩阵，其中 LAPACK 例程仅覆盖矩阵的一半。所以我使用 TRIU 将另一半归零并提取上半部分(。