调用MATLAB的内置LAPACK/BLAS例程 [英] Calling MATLAB's built-in LAPACK/BLAS routines
问题描述
我想学习如何在MATLAB中调用内置的LAPACK/BLAS例程.我有MATLAB和mex文件的经验，但实际上我不知道如何调用LAPACK或BLAS库.我已经在文件交换中找到了网关例程，因为它简化了我的通话不必为任何功能(例如此功能)编写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.
推荐答案
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函数不允许修改其输入，因此将其写为return单元格阵列中所有输入的副本(有任何修改)
我们得到:
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.实际上，文档中有一页说明
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 Files\MATLAB\R2013a\extern\lib\win64\microsoft\libmwlapack.lib"
On my 64-bit Windows, I compile the MEX-file as: 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中本征向量的符号是任意的，因此您可能会比较这两者的符号)
(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屋！
