Если вы загляните внутрь файла lapack.m из упомянутой заявки FEX, вы увидите пару примеров использования этой функции:
Пример: декомпозиция SVD с использованием DGESVD :
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 ожидают, что все аргументы будут переданы по ссылке, но MEX-функции в MATLAB не позволяют изменять их входные данные, поэтому он написан так, чтобы возвращать копии всех входных данных в массиве ячеек с любыми изменениями)
Получаем:
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
Что эквивалентно собственной функции SVD MATLAB:
[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
РЕДАКТИРОВАТЬ:
ДляДля полноты изложения ниже показан пример MEX-функции, напрямую вызывающей интерфейс Fortran подпрограммы DGESVD.
Хорошая новость заключается в том, что MATLAB предоставляет библиотеки libmwlapack и libmwblas и два соответствующих файла заголовков * 1034.* и lapack.h мы можем использовать.Фактически, в документации есть страница, объясняющая процесс вызова функций BLAS / LAPACK из MEX-файлов .
В нашем случае 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 является произвольным, поэтому вы можете получить перевернутые знаки, сравнивая их)