c-在MATLAB mexFunction中，如何将GNU科学库矩阵和向量分配给plhs

Question

我有MATLAB代码，我想称之为GNU科学库(GSL(C库，用于计算矩阵的奇异值分解。基本上，我需要实现网关MEX功能。我知道如何设置输入变量，但我该如何设置输出变量？

这是我到目前为止的代码：

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);

项目_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

Answer 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例程，则可以避免这种数据复制并提高代码的性能。