9

我想学习如何在 MATLAB 中调用内置的 LAPACK/BLAS 例程。我有 MATLAB 和 mex 文件的经验，但我实际上不知道如何调用 LAPACK 或 BLAS 库。我在文件交换中找到了简化调用的网关例程，因为我不必为任何函数（例如这个函数）编写 mex 文件。我需要任何玩具示例来学习 MATLAB 和这些内置库之间的基本消息传递。欢迎任何玩具示例，例如矩阵乘法或 LU 分解。

4

2 回答 2

10

如果您查看提到的 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-files 调用 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任意的，因此您可能会得到比较两者的翻转符号）

于 2011-06-22T13:42:11.233 回答
-2

此处解释了如何在 matlab 中使用 SVD 的一个很好/实用的示例： 将捕获的坐标转换为屏幕坐标

有关如何使用 lapack 在 Objective-C 中计算 svd 的更多信息，请参见此处计算 V from A = USVt in Objective-C with SVD from LAPACK in xcode

于 2012-08-09T21:23:07.140 回答