调用 MATLAB 的内置 LAPACK/BLAS 例程 [英] Calling MATLAB's built-in LAPACK/BLAS routines
问题描述
我想学习如何在 MATLAB 中调用内置的 LAPACK/BLAS 例程.我有使用 MATLAB 和 mex 文件的经验，但实际上我不知道如何调用 LAPACK 或 BLAS 库.我在 file exchange 中找到了网关例程，它简化了调用，因为我不必为任何函数编写 mex 文件，例如 这个.我需要任何玩具示例来学习 MATLAB 和这些内置库之间的基本消息传递.欢迎使用任何玩具示例，例如矩阵乘法或 LU 分解.
I want to learn how to call the built-in LAPACK/BLAS routines in MATLAB. I have experience in MATLAB and mex files but I've actually no idea how to call LAPACK or BLAS libraries. I've found the gateway routines in file exchange that simplifies the calls since I don't have to write a mex file for any function such as this one. I need any toy example to learn the basic messaging between MATLAB and these built-in libraries. Any toy example such as matrix multiplication or LU decomposition is welcome.
推荐答案
如果您查看 lapack.m 来自提到的 FEX 提交的文件，您将看到几个有关如何使用该功能的示例:
If you look inside the lapack.m file from the FEX submission mentioned, you will see a couple of examples on how to use the function:
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 函数:>
Which is equivalent to MATLAB's own SVD function:
[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 函数示例.
For completeness, I show below an example of a MEX-function directly calling the Fortran interface of the DGESVD routine.
好消息是 MATLAB 提供了 libmwlapack 和 libmwblas 库以及两个相应的头文件 blas.h 和 lapack.h 我们可以使用.实际上，文档中有一页解释了 从 MEX 文件调用 BLAS/LAPACK 函数.
The good news is that MATLAB provides libmwlapack and libmwblas libraries and two corresponding header files blas.h and lapack.h we can use. In fact, there is a page in the documentation explaining the process of calling BLAS/LAPACK functions from MEX-files.
在我们的例子中，lapack.h 定义了以下原型:
In our case, lapack.h defines the following prototype:
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 FilesMATLABR2013aexternlibwin64microsoftlibmwlapack.lib"
On my 64-bit Windows, I compile the MEX-file as: mex -largeArrayDims svd_lapack.c "C:Program FilesMATLABR2013aexternlibwin64microsoftlibmwlapack.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 中的特征向量的符号是任意的，因此您可能会在比较两者时得到翻转的符号)
(remember that the sign of the eigenvectors in U and V is arbitrary, so you might get flipped signs comparing the two)
这篇关于调用 MATLAB 的内置 LAPACK/BLAS 例程的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
