调用Eigen科学计算库的C++程序在Matlab环境中运行

开局一大段话（废话）

抱着提高Matlab程序运行速度的想法，打算在C++中编写耗时的计算过程，然后通过Matlab提供的Mex工具编译成mexw64程序，再用Matlab调用它。最初是打算用Matlab本身提供的科学计算库，但仔细想想，还是算了，对Matlab依赖性强不是什么好事，于是就开启了寻找开源的C++科学计算库的过程。在找库的过程中，我看到了这几个常用的库：Eigen、Ginac、Armadillo等。本来是想用Ginac的，因为它支持符号运算，但这玩意对linux系统支持得比较好，但要在windows下使用的话就需要复杂的配置安装过程，我一个学车的，真的是懒得搞这些。另外两个库都在官网上都标榜自己"快"，2选1。

当我在Eigen官网上看到这句话时，嘿嘿，直接用Eigen了，在windows下使用Eigen库，只需要复制粘贴，再#include就完了，实在是方便。

Eigen doesn’t have any dependencies other than the C++ standard library.If you just want to use Eigen, you can use the header files right away. There is no binary library to link to, and no configured header file. Eigen is a pure template library defined in the headers.

然而，今天发现用C++写的程序运行速度还比不上Matlab自身的运行速度。只有当处理的矩阵尺寸比较小时，C++写的程序才能在运行速度上快人一步。然并卵噢。
测试结果：

>> timeSpendCplus(2^9,1000);
矩阵行数和列数分别为512,512

C++程序1000个矩阵相乘所花费时间为25.397602

>> timeSpendMat(2^9,1000);
矩阵行数和列数分别为512,512

m程序1000个矩阵相乘所花费时间为3.279014

1 windows下使用VScode调用Eigen库

• 第一步，官网下载Eigen库，然后解压，比方说解压到C:\usr\peter\Desktop\下，得到文件夹eigen-3.x.x.

• 第二步，使用的时候包含头文件 #include <Eigen\Dense>

• 第三步，编译时使用语句：g++ demo.cpp -o demo -IC:\usr\peter\Desktop\eigen-3.x.x\

• 第四步，使用./demo运行第三步生成的demo.exe，查看结果是否符合预期

1.1 实例

demo.cpp:

//demo.cpp
#include <iostream>
#include <Eigen/Dense>
 
using Eigen::MatrixXd;
 
int main()
{
  MatrixXd m(2,2);
  m(0,0) = 3;
  m(1,0) = 2.5;
  m(0,1) = -1;
  m(1,1) = m(1,0) + m(0,1);
  std::cout << m << std::endl;
}

编译：

g++ demo.cpp -o demo -IC:\usr\peter\Desktop\eigen-3.x.x\

2 mexFunction中使用Eigen库

mexFunction是Matlab使用mex编译.cpp和.c文件的入口函数，就相当于主函数。

2.1 mexFunction使用实例

demo.cpp：

//demo.cpp

#include "mex.h"//包含调用matlab函数需要的头文件

// lhs = left hand side
// rhs = right hand side
void mexFunction(int nlhs, mxArray * plhs[], int nrhs, const mxArray *prhs[])
{
    //prhs[0]指向第一个输入参数
    //prhs[1]指向第二个输入参数
    //每个输入参数的都是mxArray类型
    //mxArray数组里的每个元素都是double数据类型
    
  
    //假设只允许输入一个参数
    if(nrhs != 1){ mexErrMsgTxt("please input only one argument");}
    //mexErrorMsgTxt是Matlab里的函数，真是又长又臭
    //注意使用mexErrorMsgTxt程序就结束了，且不会自动回收内存，所以该函数应该放在回收内存之后
    //不然会造成内存泄漏
    
    double *x;//定义一个指向数组里的double元素的指针
    x = mxGetPr(prhs[0]);//mxGetPr返回数组第一个元素的指针
    
    double *y;
    //假设该函数只有一个返回参数
    //将返回指针指向一个double矩阵
    int rows=1,cols=1;//行列数均为3
    plhs[0] = mxCreateDoubleMatrix(rows,cols,mxREAL);
    y = mxGetPr(plhs[0]);
    
    //以上获得了来自Matlab的输入数据并创建了返回给Matlab的数据结构
    
    //计算部分：
    //假设mexFunction要实现的功能是output = 2* input
    
    y[0]=2*x[0];//OK
    
}

输入以下语句编译得到demo.mexw64文件：

mex demo.cpp  -output demo

测试：

>> test(80)

ans =

   160

另外Matlab传过来的矩阵是按列排列的，比如一个3x3矩阵的第2行第一个元素，传至C++后的索引号是[1]，第1行第2列元素的索引号是[3]。

2.2 使用Eigen库的mexFunction函数

demo.cpp：写得乱，看着玩吧

//demo.cpp

#include "mex.h"
#include <Eigen\Dense>

using namespace Eigen;



//因为两个矩阵不能直接赋值，写函数使两个结构互相转换
void matToEigen(MatrixXd & m,const mxArray *prhs, int mxRows,int mxCols);
void EigenTomat(double *y ,const MatrixXd & m);//只在最后传递返回值调用

//Matrixd类是Eigen库的类
//mx和mex开头的函数都是Matlab的函数

void mexFunction(int nlhs, mxArray * plhs[], int nrhs, const mxArray *prhs[])
{
    
    //获得matlab矩阵的函数和列数
    int mxRows,mxCols;
    mxRows = mxGetM(prhs[0]);
    mxCols = mxGetN(prhs[0]);
    
    //Eigen库对象
    MatrixXd m0(mxRows,mxCols);
    MatrixXd & m1 =m0;
    MatrixXd m2(mxRows,mxCols);
    MatrixXd & m3 = m2;
  
    matToEigen(m1,prhs[0],mxRows,mxCols);//m1接收到了matlab传来的矩阵
    
    // 使用Eigen库进行计算
    m3 = 2*m0;
    
   	//创建一个matlab接收用的数组
    double *y;
    plhs[0] = mxCreateDoubleMatrix(mxRows,mxCols,mxREAL);
    y = mxGetPr(plhs[0]);
    
    EigenTomat(y ,m3);//计算结果发回matlab
    
}


void matToEigen(MatrixXd & m,const mxArray *prhs_, int mxRows,int mxCols){
    double *x;

    
    x = mxGetPr(prhs_);//指向数组元素
    
    //赋值给Eigen矩阵
   	for(int k=0; k<mxRows; k++){
        for(int l=0; l<mxCols; l++)
        {
            m(k,l) = x[mxRows*l+k];
        }
    }
}



void EigenTomat(double *y ,const MatrixXd & m){
    int rows,cols,index;
    rows = m.rows();
    cols = m.cols(); 

    for(int r=0; r<rows; r++){
        for(int c=0; c<cols; c++)
        {
            y[rows*c+r] = m(r,c);//赋值给Y,传到Matlab中
        }
    }
    
}

输入以下语句编译得到demo.mexw64文件：

mex demo.cpp -output demo -IC:\usr\peter\Desktop\eigen-3.x.x\

测试：

>> demo([2,3;0,1])

ans =

     4     6
     0     2

另外Matlab传过来的矩阵是按列排列的，比如一个3x3矩阵的第2行第一个元素，传至C++后的索引号是[1]，第1行第2列元素的索引号是[3]。

3 可能会遇到的小坑

• 入口函数书写错误

  void mexFunction(int nlhs, mxArray * plhs[], int nrhs, const mxArray *prhs[])
  

  //如果第四个参数遗漏了const

  //则会报错：
  错误使用 mex
  Cannot export mexFunction: symbol not defined
  collect2.exe: error: ld returned 1 exit status

• 忘记使用mxCreateDoubleMatrix()或相似类型声明传送到Matlab的数据类型，编译没有问题，但到matlab中运行会报错：