我正在用 mex(使用 C)重新编写一段 MATLAB 代码。到目前为止,我的 C 版本的 MATLAB 代码大约是 MATLAB 代码的两倍。现在我有三个问题,都与下面的代码有关:
- 我怎样才能加快这段代码的速度?
- 你发现这段代码有什么问题吗?我问这个是因为我不太了解 mex,而且我也不是 C 大师 ;-) ...我知道应该在代码中进行一些检查(例如,如果在使用时仍有堆空间
realloc
,但为了简单起见,我暂时将其保留) - 有没有可能,MATLAB 优化得如此之好,以至于我真的无法获得比 C 语言快两倍以上的代码...?
代码应该或多或少与平台无关(Win、Linux、Unix、Mac、不同的硬件),所以我不想使用汇编程序或特定的线性代数库。所以这就是为什么我自己对员工进行编程...
#include <mex.h>
#include <math.h>
#include <matrix.h>
void mexFunction(
int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
{
double epsilon = ((double)(mxGetScalar(prhs[0])));
int strengthDim = ((int)(mxGetScalar(prhs[1])));
int lenPartMat = ((int)(mxGetScalar(prhs[2])));
int numParts = ((int)(mxGetScalar(prhs[3])));
double *partMat = mxGetPr(prhs[4]);
const mxArray* verletListCells = prhs[5];
mxArray *verletList;
double *pseSum = (double *) malloc(numParts * sizeof(double));
for(int i = 0; i < numParts; i++) pseSum[i] = 0.0;
float *tempVar = NULL;
for(int i = 0; i < numParts; i++)
{
verletList = mxGetCell(verletListCells,i);
int numberVerlet = mxGetM(verletList);
tempVar = (float *) realloc(tempVar, numberVerlet * sizeof(float) * 2);
for(int a = 0; a < numberVerlet; a++)
{
tempVar[a*2] = partMat[((int) (*(mxGetPr(verletList) + a))) - 1] - partMat[i];
tempVar[a*2 + 1] = partMat[((int) (*(mxGetPr(verletList) + a))) - 1 + lenPartMat] - partMat[i + lenPartMat];
tempVar[a*2] = pow(tempVar[a*2],2);
tempVar[a*2 + 1] = pow(tempVar[a*2 + 1],2);
tempVar[a*2] = tempVar[a*2] + tempVar[a*2 + 1];
tempVar[a*2] = sqrt(tempVar[a*2]);
tempVar[a*2] = 4.0/(pow(epsilon,2) * M_PI) * exp(-(pow((tempVar[a*2]/epsilon),2)));
pseSum[i] = pseSum[i] + ((partMat[((int) (*(mxGetPr(verletList) + a))) - 1 + 2*lenPartMat] - partMat[i + (2 * lenPartMat)]) * tempVar[a*2]);
}
}
plhs[0] = mxCreateDoubleMatrix(numParts,1,mxREAL);
for(int a = 0; a < numParts; a++)
{
*(mxGetPr(plhs[0]) + a) = pseSum[a];
}
free(tempVar);
free(pseSum);
}
所以这是改进版,比MATLAB版快12倍左右。转换的事情仍然占用了很多时间,但我暂时放弃了,因为我必须为此在 MATLAB 中更改一些东西。所以首先关注剩下的C代码。您是否在以下代码中看到更多潜力?
#include <mex.h>
#include <math.h>
#include <matrix.h>
void mexFunction(
int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
{
double epsilon = ((double)(mxGetScalar(prhs[0])));
int strengthDim = ((int)(mxGetScalar(prhs[1])));
int lenPartMat = ((int)(mxGetScalar(prhs[2])));
double *partMat = mxGetPr(prhs[3]);
const mxArray* verletListCells = prhs[4];
int numParts = mxGetM(verletListCells);
mxArray *verletList;
plhs[0] = mxCreateDoubleMatrix(numParts,1,mxREAL);
double *pseSum = mxGetPr(plhs[0]);
double epsilonSquared = epsilon*epsilon;
double preConst = 4.0/((epsilonSquared) * M_PI);
int numberVerlet = 0;
double tempVar[2];
for(int i = 0; i < numParts; i++)
{
verletList = mxGetCell(verletListCells,i);
double *verletListPtr = mxGetPr(verletList);
numberVerlet = mxGetM(verletList);
for(int a = 0; a < numberVerlet; a++)
{
int adress = ((int) (*(verletListPtr + a))) - 1;
tempVar[0] = partMat[adress] - partMat[i];
tempVar[1] = partMat[adress + lenPartMat] - partMat[i + lenPartMat];
tempVar[0] = tempVar[0]*tempVar[0] + tempVar[1]*tempVar[1];
tempVar[0] = preConst * exp(-(tempVar[0]/epsilonSquared));
pseSum[i] += ((partMat[adress + 2*lenPartMat] - partMat[i + (2*lenPartMat)]* tempVar[0]);
}
}
}
最佳答案
您不需要分配 pseSum 供本地使用,然后再将数据复制到输出。您可以简单地分配一个 MATLAB 对象并获取指向内存的指针:
plhs[0] = mxCreateDoubleMatrix(numParts,1,mxREAL); pseSum = mxGetPr(plhs[0]);
因此您不必将 pseSum 初始化为 0,因为 MATLAB 已经在 mxCreateDoubleMatrix 中完成了它。
将内循环中的mxGetPr全部去掉,赋值给之前的变量。
考虑在 MATLAB 中使用 int32 或 uint32 数组,而不是将 double 转换为整数。将 double 转换为 int 是昂贵的。内部循环计算看起来像
tempVar[a*2] = partMat[somevar[a] - 1] - partMat[i];
你在你的代码中使用了这样的结构
((int) (*(mxGetPr(verletList) + a)))
您这样做是因为 varletList 是一个“ double ”数组(在 MATLAB 中默认情况下就是这种情况),它包含整数值。相反,您应该使用整数数组。在 MATLAB 中调用 mex 文件类型之前:
varletList = int32(varletList);
那么你就不需要上面的类型转换为 int 了。你会简单地写
((int*)mxGetData(verletList))[a]
或者更好的是,更早分配
somevar = (int*)mxGetData(verletList);
以后再写
somevar[a]
在所有循环之前预计算 4.0/(pow(epsilon,2) * M_PI)!这是一个昂贵的常量。
pow((tempVar[a*2]/epsilon),2)) 就是 tempVar[a*2]^2/epsilon^2。您之前计算 sqrt(tempVar[a*2]) 。你为什么现在平方它?
一般不使用 pow(x, 2)。就写x*x
我会在参数上添加一些合理性检查,尤其是当您需要整数时。要么使用 MATLAB 的 int32/uint32 类型,要么检查你实际得到的是一个整数。
在新代码中编辑
在循环之前计算 -1/epsilonSquared 并计算 exp(minvepssq*tempVar[0])。请注意,结果可能略有不同。取决于您的需要,但如果您不关心具体的操作顺序,那就去做吧。
定义一个寄存器变量 preSum_r 并使用它对内部循环中的结果求和。循环后将其分配给 preSum[i]。如果你想要更多乐趣,可以使用 SSE streaming store(_mm_stream_pd compiler intrinsic)将结果写入内存。
移除 double 到 int 的转换
很可能无关紧要,但请尝试将 tempVar[0/1] 更改为普通变量。无关紧要,因为编译器应该为你做那件事。但同样,这里不需要数组。
使用 OpenMP 并行化外部循环。简单(至少是不考虑 NUMA 架构的数据布局的最简单版本),因为迭代之间没有依赖性。
关于c - 如何加速这个墨西哥代码?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/12272218/