在 MATLAB mexFunction 中，如何将 GNU 科学库矩阵和向量分配给 plhs？
Posted
技术标签:
【中文标题】在 MATLAB mexFunction 中，如何将 GNU 科学库矩阵和向量分配给 plhs？【英文标题】：In the MATLAB mexFunction how do you assign GNU Scientific Library matrices and vectors to plhs? 【发布时间】：2021-03-17 15:46:02 【问题描述】：我有 MATLAB 代码，我想调用 GNU 科学库 (GSL) C 库来计算矩阵的 singular value decomposition。基本上，我需要实现gateway MEX function。我知道如何设置输入变量，但是如何设置输出变量？
到目前为止，这是我的代码：
build.m
mex -v -L/home/miran045/reine097/c-libs/lib/ -R2017b project_m.c -lgsl -lgslcblas -lm
call_c_from_matlab.m
A_data = [1.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0];
[U, S, V] = project_m(4, 5, A_data);
project_m_from_c
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#include "mex.h"
void run_svd(const size_t M, const size_t N, double A_data[], gsl_matrix * B, gsl_matrix * V, gsl_vector * S, gsl_vector * work)
gsl_linalg_SV_decomp(B, V, S, work);
static double
get_double(const mxArray *array_ptr)
double *pr;
pr = mxGetPr(array_ptr);
double result = (double) *pr;
return result;
void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )
int i, j;
/* You should put code here to check nrhs, nlhs, prhs[0], prhs[1], prhs[2] */
int m = (int) get_double(prhs[0]);
int n = (int) get_double(prhs[1]);
double *mexArray = mxGetPr(prhs[2]);
int Bm, Bn; /* added */
gsl_matrix * B;
gsl_matrix * V;
gsl_vector * S;
gsl_vector * work;
gsl_matrix_view A = gsl_matrix_view_array(mexArray, m, n);
if (n > m)
Bm = n; Bn = m; /* added */
work = gsl_vector_alloc(m);
gsl_matrix_transpose_memcpy(B, &A.matrix);
else
Bm = m; Bn = n; /* added */
work = gsl_vector_alloc(n);
gsl_matrix_memcpy(B, &A.matrix);
plhs[0] = mxCreateDoubleMatrix( Bm, Bn, mxREAL );
plhs[1] = mxCreateDoubleMatrix( Bn, 1, mxREAL );
plhs[2] = mxCreateDoubleMatrix( Bn, Bn, mxREAL );
run_svd(m, n, mexArray, plhs[0], plhs[2], plhs[1], work);
gsl_vector_free(work);
但是我收到以下错误：
gsl: ../gsl/gsl_matrix_double.h:278: ERROR: first index out of range
Default GSL error handler invoked.
matlab&
[2] 130182
[1] Killed matlab
reine097@native [/panfs/roc/groups/4/miran045/reine097] % MATLAB is selecting SOFTWARE OPENGL rendering.
MATLAB: smart_ptr.cpp:82: matrix::unique_mxarray_ptr matrix::from_matlab(mxArray* const&): Assertion `Input to from_matlab must have temporary or local scope' failed.
matlab&
[3] 115869
[2] Segmentation fault (core dumped) matlab
reine097@native [/panfs/roc/groups/4/miran045/reine097] % MATLAB is selecting SOFTWARE OPENGL rendering.
reine097@native [/panfs/roc/groups/4/miran045/reine097] % gsl: svd.c:79: ERROR: length of workspace must match second dimension of matrix A
Default GSL error handler invoked.
matlab&
[4] 89615
[3] Killed matlab
【问题讨论】：
【参考方案1】：您的代码存在多个问题：
您没有检查输入是否完全符合预期（实数、双精度、非稀疏等） 您正在将未分配 MATLAB API 函数的内存附加到 mxArray 在将指针附加到 mxArray 后立即释放指针后面的内存 在使用 plhs[ ] mxArray 变量之前不要创建它们通常，在使用 mex 例程时，会按顺序执行以下操作：
检查输入是否完全符合预期 创建输出 plhs[ ] mxArray 变量 获取指向 plhs[ ] mxArray 变量的数据区的指针 将这些指针传递给计算例程在您的情况下，您可能出于某种原因想要使用 gsl_matrix_alloc( ) 例程进行内存分配，而不是使用 MATLAB API 内存分配函数。很好，但是如果你这样做，那么你必须将结果复制到 plhs[ ] mxArray 数据区。您不能直接附加指针，因为这会搞砸 MATLAB 内存管理器并导致下游崩溃。但在所有情况下，您都需要先创建 plhs[ ] mxArray 变量。例如，如果我们只是使用您当前的内存分配方案，您将需要复制数据：
编辑 gsl 数据复制
/* Create MATLAB mxArray from gsl_vector */
mxArray *gslvector2MATLAB(gsl_vector *V)
mxArray *mx;
double *data, *target;
size_t n, stride;
if( V )
n = V->size;
data = V->data;
stride = V->stride;
mx = mxCreateDoubleMatrix(n,1,mxREAL);
target = (double *) mxGetData(mx);
while( n-- )
*target++ = *data;
data += stride;
else
mx = mxCreateDoubleMatrix(0,0,mxREAL);
return mx;
/* Create MATLAB mxArray from gsl_matrix */
mxArray *gslmatrix2MATLAB(gsl_matrix *M)
mxArray *mx;
double *data, *target;
size_t i, j, m, n;
if( M )
m = M->size1;
n = M->size2;
data = M->data;
mx = mxCreateDoubleMatrix(m,n,mxREAL);
target = (double *) mxGetData(mx);
for( i=0; i<m; i++ )
for( j=0; j<n; j++ )
target[i+j*m] = data[j]; /* tranpose */
data += M->tda;
else
mx = mxCreateDoubleMatrix(0,0,mxREAL);
return mx;
void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )
int i;
/* You should put code here to check nrhs, nlhs, prhs[0], prhs[1], prhs[2] */
int m = (int) get_double(prhs[0]);
int n = (int) get_double(prhs[1]);
double *mexArray = mxGetPr(prhs[2]);
gsl_matrix * B;
gsl_matrix * V;
gsl_vector * S;
gsl_vector * work;
gsl_matrix_view A = gsl_matrix_view_array(mexArray, m, n);
if (n > m)
B = gsl_matrix_alloc(n, m);
V = gsl_matrix_alloc(m, m);
S = gsl_vector_alloc(m);
work = gsl_vector_alloc(m);
gsl_matrix_transpose_memcpy(B, &A.matrix);
else
B = gsl_matrix_alloc(m, n);
V = gsl_matrix_alloc(n, n);
S = gsl_vector_alloc(n);
work = gsl_vector_alloc(n);
gsl_matrix_memcpy(B, &A.matrix);
run_svd(m, n, mexArray, B, V, S, work);
nlhs = 3; // delete this line
mxSetPr(plhs[0], B); // delete this line
mxSetPr(plhs[1], S); // delete this line
mxSetPr(plhs[2], V); // delete this line
/* call custom routines to copy data as transpose */
plhs[0] = gslmatrix2MATLAB(B);
plhs[1] = gslvector2MATLAB(S);
plhs[2] = gslmatrix2MATLAB(V);
gsl_matrix_free(B);
gsl_matrix_free(V);
gsl_vector_free(S);
gsl_vector_free(work);
我不熟悉 gsl_matrix 或 gsl_vector 类型，因此无法建议您如何复制上述数据。通常，您可以从 gsl 类型中获取数据指针，然后 memcpy 数据。如果您遵循我上面概述的步骤并首先创建 plhs[ ] 变量，然后将它们的数据指针直接传递给您的 svd 例程，您可以避免这种数据复制并提高代码的性能。
【讨论】：
关于这一点：/* 在此处添加代码以从 B 复制到 mxGetData(plhs[0]) */ 我知道如何从 B 获取第 i 行和第 j 列。你这样做：gsl_matrix_get(B, i, j)。但是如何设置mxGetData(plhs[0])的第i行第j列呢？ 我查看了 gsl 向量和矩阵类型的文档，它看起来像一个简单的结构，其中包含一个名为“data”的字段，它是一个指向 double 的指针。所以看起来你可以这样做： memcpy(mxGetData(plhs[0]),B->data,Bm*Bn); S 和 V 也是如此。 好吧，我进一步阅读了 gsl 文档，发现矩阵以行优先顺序存储（如 C）。但是 MATLAB 以列优先顺序存储矩阵（如 Fortran）。所以上面的 memcpy() 不会直接起作用。您基本上需要在 C 代码中或在 MATLAB 端获得矩阵后转置矩阵。如果您在这部分需要帮助，请告诉我。 请注意，如果 tda 与 gsl 矩阵中的 size2 不匹配，则基础数据将不连续，您必须逐行复制循环中的数据。 这部分我需要一些帮助，如何设置 mxGetData(plhs[0]) 的第 i 行和第 j 列元素？以上是关于在 MATLAB mexFunction 中，如何将 GNU 科学库矩阵和向量分配给 plhs？的主要内容，如果未能解决你的问题，请参考以下文章
我应该如何将数据传递给 mexFunction 的 mxArray *plhs[] 以便在 Matlab 中获得其输出？
在 Visual Studio 中构建 MATLAB mex 文件会给出“函数 mexFunction 中引用的 LNK2019 未解析的外部符号 _mexPrintf”？