在 MATLAB mexFunction 中，如何将 GNU Scientific Library 矩阵和向量分配给 plhs？

问题描述

我有 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,2.0,3.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,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,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

解决方法

您的代码存在多个问题：

• 您没有检查输入是否完全符合预期（真实、双重、非稀疏等）
• 您将在没有 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,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,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,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,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 例程，则可以避免这种数据复制并提高代码的性能。