最近在超算上做大规模的蒙特卡洛仿真，发现自己的matlab代码非常耗时，利用matlab的profiler工具分析代码，发现耗时的原因是一个没办法vectorizing的函数海量执行for循环。查找matlab优化代码方法，网上给的建议是关键代码用c语言实现，然后混合编程。网上搜索混合编程方法，发现大多数讲的都太偏mex的使用了，还有人说mex中的变量类型的（如mwSize），搞的我以为标准c语言不能再用了。其实我已经具有该关键函数的c语言实现，需要的只是一个接口函数而已。

代码如下：

在一个c文件中，保留原来的所有代码，并加入#include "mex.h", 以及接口函数

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

其中接口函数的返回类型和输入参数都是固定的。在接口函数中加入自己的代码，规定好输入输出，中间就可以调用自己已有的函数了。改好之后，用mex 编译该文件，文件名就是matlab中调用的函数名。例如，我的c文件名为：csearch.c，改写完之后，在matlab命令行中输入“mex csearch.c”，如果编译成功，就可以在matlab中使用csearch函数了。输入输出的规定见接口函数的代码注释。

#include "mex.h"

#include "stdlib.h"

#include "stdio.h"

#include "math.h"

#define LOOPMAX     10000           /* maximum count of search loop */

#define SGN(x)      ((x)<=0.0?-1.0:1.0)

#define ROUND(x)    (floor((x)+0.5))

#define SWAP(x,y)   do {double tmp_; tmp_=x; x=y; y=tmp_;} while (0)

/* new matrix ------------------------------------------------------------------

* allocate memory of matrix 

* args   : int    n,m       I   number of rows and columns of matrix

* return : matrix pointer (if n<=0 or m<=0, return NULL)

*-----------------------------------------------------------------------------*/

double *mat(int n, int m)

{

    double *p;

    ... 

    ...

    return p;

}

/* zero matrix -----------------------------------------------------------------

* generate new zero matrix

* args   : int    n,m       I   number of rows and columns of matrix

* return : matrix pointer (if n<=0 or m<=0, return NULL)

*-----------------------------------------------------------------------------*/

double *zeros(int n, int m)

{

    double *p;

    ...

    ...

    return p;

}

int search(int n, int m, const double *L, const double *D,

                  const double *zs, double *zn, double *s)

{

    int i,j,k,c,nn=0,imax=0;

    double newdist,maxdist=1E99,y;

    double *S=zeros(n,n),*dist=mat(n,1),*zb=mat(n,1),*z=mat(n,1),*step=mat(n,1);

       ...

       ...

       ...

free(S); free(dist); free(zb); free(z); free(step);

    if (c>=LOOPMAX) {

        fprintf(stderr,"%s : search loop count overflow\n",__FILE__);

        return -1;

    }

    return 0;

}

//接口函数：nlhs--输出参数个数；*plhs[]--输出参数组；nrhs--输入参数个数；*prhs[]--输入参数组

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

{

// [afixed,sqnorm]=csearch(ahat,L,D,ncands) %matlab函数。在matlab中，我会这样调用这个函数，所以输入参数有4个：ahat--nx1维列向量；L--nxn维矩阵；D--nx1维列向量; ncands--int 类型数值

// 输出参数：afixed--nxm维矩阵；sqnorm--mx1维列向量

  double *L,*D,*zs,*s,*zn; //定义指向输入输出的指针

  int n,m=2;

 //得到输入参数

  n = mxGetM(prhs[0]);// 得到第一个输入参数ahat的行数

  zs= mxGetPr(prhs[0]);//得到指向第一个输入参数ahat的指针

  L = mxGetPr(prhs[1]);//得到指向第二个输入参数L的指针

  D = mxGetPr(prhs[2]);//得到指向第三个输入参数D的指针

  m = mxGetScalar(prhs[3]);//得到指向第四个输入参数ncands的数值大小

//定义输出参数 [afixed,sqnorm]=csearch(ahat,L,D,ncands) 

  plhs[0] = mxCreateDoubleMatrix(n,m, mxREAL); // 创建输出参数1：afixed--nxm维矩阵

  plhs[1] = mxCreateDoubleMatrix(1,m, mxREAL); // 创建输出参数2：sqnorm--mx1维列向量

 //将指针zn、s分别指向输出参数1、参数2的地址，以便将参数1、参数2代入到c函数search中计算。

  zn=(double*)mxGetData(plhs[0]);

  s=(double*)mxGetData(plhs[1]);

//将参数代入到search函数中计算。其中zn和s为输出参数。

  search(n, m, L, D, zs, zn,s);//c语言函数。

//计算完成后，输出参数plhs[0] 、plhs[1] 的地址中已经具有了数值。在matlab中调用csearch将会得到期望的

//输出。

}

#######################################################################

最后，展示使用混合编程带来的效率提升：左边是使用了混合编程，matlab+c语言，右边是没有进行混合编程，即纯matlab编程。

运行了100次，matlab+c用了0.114秒，纯matlab用了10.886s,。

如果运行1E6次，用时将分别为： 1140s=19分钟  vs  108860s=30.24 小时