调用MATLAB的内置LAPACK/BLAS例程

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.039753s = 2.1354 0.88509 0.27922V = -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 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.0199S = 1.0337 0.5136 0.0811VT = -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.0199ans = 1.0337 0.5136 0.0811ans = -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例程的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！ 上岸，阿里云！