调用MATLAB的内置LAPACK/BLAS例程

Question

我想学习如何在MATLAB中调用内置的LAPACK/BLAS例程.我有MATLAB和mex文件的经验但我实际上不知道如何调用LAPACK或BLAS库.我发现在网关程序文件交换,简化了电话,因为我没有写任何功能的MEX文件,如这一个.我需要任何玩具示例来学习MATLAB和这些内置库之间的基本消息传递.任何玩具示例,如矩阵乘法或LU分解都是受欢迎的.

Answer 1

如果你从提到的FEX提交中查看lapack.m文件,你会看到几个关于如何使用该函数的例子:

示例:使用DGESVD进行 SVD分解:

X = rand(4,3);
[m,n] = size(X);
C = lapack('dgesvd', ...
     'A', 'A', ...           % compute ALL left/right singular vectors
      m, n, X, m, ...        % input MxN matrix
      zeros(n,1), ...        % output S array
      zeros(m), m, ...       % output U matrix
      zeros(n), n, ....      % output VT matrix
      zeros(5*m,1), 5*m, ... % workspace array
      0 ...                  % return value
);
[s,U,VT] = C{[7,8,10]};      % extract outputs
V = VT';

(注意:我们将这些虚拟变量用于输出变量的原因是因为Fortran函数希望所有参数都通过引用传递,但MATLAB中的MEX函数不允许修改它们的输入,因此它被编写为返回所有输入的副本任何修改的单元格数组)

我们得到:

U =
     -0.44459      -0.6264     -0.54243       0.3402
     -0.61505     0.035348      0.69537      0.37004
     -0.41561     -0.26532      0.10543     -0.86357
     -0.50132      0.73211     -0.45948    -0.039753
s =
       2.1354
      0.88509
      0.27922
V =
     -0.58777      0.20822     -0.78178
      -0.6026     -0.75743      0.25133
     -0.53981      0.61882      0.57067

这相当于MATLAB自己的SVD功能:

[U,S,V] = svd(X);
s = diag(S);

这给了:

U =
     -0.44459      -0.6264     -0.54243       0.3402
     -0.61505     0.035348      0.69537      0.37004
     -0.41561     -0.26532      0.10543     -0.86357
     -0.50132      0.73211     -0.45948    -0.039753
s =
       2.1354
      0.88509
      0.27922
V =
     -0.58777      0.20822     -0.78178
      -0.6026     -0.75743      0.25133
     -0.53981      0.61882      0.57067

编辑:

为了完整起见,我在下面显示了一个直接调用DGESVD例程的Fortran接口的MEX函数示例.

好消息是,MATLAB提供libmwlapack与libmwblas图书馆和两个相应的头文件blas.h,并lapack.h可以使用.实际上,文档中有一个页面解释了从MEX文件调用BLAS/LAPACK函数的过程.

在我们的例子中,lapack.h定义以下原型:

extern void dgesvd(char *jobu, char *jobvt, 
  ptrdiff_t *m, ptrdiff_t *n, double *a, ptrdiff_t *lda,
  double *s, double *u, ptrdiff_t *ldu, double *vt, ptrdiff_t *ldvt,
  double *work, ptrdiff_t *lwork, ptrdiff_t *info);

svd_lapack.c

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

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    mwSignedIndex m, n, lwork, info=0;
    double *A, *U, *S, *VT, *work;
    double workopt = 0;
    mxArray *in;

    /* verify input/output arguments */
    if (nrhs != 1) {
        mexErrMsgTxt("One input argument required.");
    }
    if (nlhs > 3) {
        mexErrMsgTxt("Too many output arguments.");
    }
    if (!mxIsDouble(prhs[0]) || mxIsComplex(prhs[0])) {
        mexErrMsgTxt("Input matrix must be real double matrix.");
    }

    /* duplicate input matrix (since its contents will be overwritten) */
    in = mxDuplicateArray(prhs[0]);

    /* dimensions of input matrix */
    m = mxGetM(in);
    n = mxGetN(in);

    /* create output matrices */
    plhs[0] = mxCreateDoubleMatrix(m, m, mxREAL);
    plhs[1] = mxCreateDoubleMatrix((m<n)?m:n, 1, mxREAL);
    plhs[2] = mxCreateDoubleMatrix(n, n, mxREAL);

    /* get pointers to data */
    A = mxGetPr(in);
    U = mxGetPr(plhs[0]);
    S = mxGetPr(plhs[1]);
    VT = mxGetPr(plhs[2]);

    /* query and allocate the optimal workspace size */
    lwork = -1;
    dgesvd("A", "A", &m, &n, A, &m, S, U, &m, VT, &n, &workopt, &lwork, &info);
    lwork = (mwSignedIndex) workopt;
    work = (double *) mxMalloc(lwork * sizeof(double));

    /* perform SVD decomposition */
    dgesvd("A", "A", &m, &n, A, &m, S, U, &m, VT, &n, work, &lwork, &info);

    /* cleanup */
    mxFree(work);
    mxDestroyArray(in);

    /* check if call was successful */
    if (info < 0) {
        mexErrMsgTxt("Illegal values in arguments.");
    } else if (info > 0) {
        mexErrMsgTxt("Failed to converge.");
    }
}

在我的64位Windows上,我将MEX文件编译为: mex -largeArrayDims svd_lapack.c "C:\Program Files\MATLAB\R2013a\extern\lib\win64\microsoft\libmwlapack.lib"

这是一个测试:

>> X = rand(4,3);
>> [U,S,VT] = svd_lapack(X)
U =
   -0.5964    0.4049    0.6870   -0.0916
   -0.3635    0.3157   -0.3975    0.7811
   -0.3514    0.3645   -0.6022   -0.6173
   -0.6234   -0.7769   -0.0861   -0.0199
S =
    1.0337
    0.5136
    0.0811
VT =
   -0.6065   -0.5151   -0.6057
    0.0192    0.7521   -0.6588
   -0.7949    0.4112    0.4462

与

>> [U,S,V] = svd(X);
>> U, diag(S), V'
U =
   -0.5964    0.4049    0.6870    0.0916
   -0.3635    0.3157   -0.3975   -0.7811
   -0.3514    0.3645   -0.6022    0.6173
   -0.6234   -0.7769   -0.0861    0.0199
ans =
    1.0337
    0.5136
    0.0811
ans =
   -0.6065   -0.5151   -0.6057
    0.0192    0.7521   -0.6588
   -0.7949    0.4112    0.4462

(记住,特征向量的符号是U和V任意的,所以你可能会得到翻转的符号来比较两者)