调用MATLAB内置的LAPACK/BLAS例程

Question

调用MATLAB内置的LAPACK/BLAS例程

9

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

- petrichor

1

为什么要使用FORTRAN？我不确定，但我一直认为LAPACK/BLAS已经在30年前移植到了C/C++上，并且MATLAB使用了这个实现！ - Mikhail Poda

嗯，你看我对此很陌生 :( 当我查找文档时，我看到了FORTRAN代码。让我相应地更新问题。 - petrichor

3

@Mikhail, @Ismail：这些库本身是用Fortran编写的，调用约定也是Fortran的。你们说得对，有一些C/C++的封装器，其中最简单的是CBLAS/CLAPACK，还有一堆用于C++的模板/通用/其他东西。 - ev-br

2个回答

-2

一个良好/实用的Matlab中使用SVD的例子可以在此处找到：将捕获的坐标转换为屏幕坐标更多有关如何在Objective-C中使用LAPACK计算SVD的内容可以在此处找到：使用来自Xcode的SVD从LAPACK中计算A=USVt的V

- fellowworldcitizen

尊敬的fellowworldcitizen，欢迎来到SO。该问题询问如何从MATLAB调用内置的LAPACK例程。您的答案提供了很好的参考资料，但并不是该问题的真正答案。 - petrichor

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Amro · Accepted Answer

如果您查看来自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

vs.

：意为“对比”或“与”，用于表示两个或多个事物之间的比较或对比。

>> [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中特征向量的符号是任意的，因此在比较两者时可能会出现翻转的符号)