Matlab 和 C++ （2）
Matlab中用到c和汇编的一个例子：老板让在非线性拟和时要计算一个评价函数：N= sum( diag( B'*B ) ) ; B= (A*S) .* log( (A*S) ./ X ) - A*S + X ; 因为A是一个很大的矩阵(1024*1024),所以，这个评价函数计算很慢。最慢的一步还是log。所以，把整个评价函数用c重新写了。
其中的log用的是汇编，发现还比较有用，就单独做了一个myLOG.m以代替matlab的log（计算自然对数ln值）
检验：
>> a=rand(1000,1000,10)*10;
>> tic;b=log(a);toc
Elapsed time is 2.521904 seconds.
>> tic;mylog(a);toc （为了节约内存，我们直接将输入的矩阵作为输出）
Elapsed time is 1.531487 seconds.
>> sqrt(sum(sum(sum((a-b).^2)))/1000/1000/10) （误差）
ans =
1.2344e-016
误差的来源主要是当计算的数值比较小的时候（1e-3），我们的汇编计算没有将计算值先修正为2附近的值。但是如果计算的值不是十分小，精度还是很好的。
主要的进步还是减少了时间。
不过，matlab在单步运算上还是比较优化的。假如不是汇编计算，用真的c算log，速度和matlab差不多。（文件中间定义了#define myln_sf(__a) log(__a)，可以用来比较c的log运算）。
另外，这个汇编在计算时，没有对数据检验，所以，如果用了负数，计算的速度慢不算，结果还不对，所以用的时候要注意。
用同样的方法还可以优化exp() sin() cos()，这些函数都很有用。尤其在计算fft的时候。
以后会贴一个完整的myfft，为的是独立于matlab。
//-----------------------------myLOG.m-------------------------------------------
#include <stdio.h>
//#include <time.h>
//#include <sys\timeb.h>
#include "mex.h"
#include "math.h"
//---------------------------------------------------------------------------
static long double dM_LOG2E = 1.4426950408889634073599246810019;
static long double dM_LN2 = 0.69314718055994530941723212145818;
static long double dONE = 1;
//---------------------------------------------------------------------------
//#define myln_sf(__a) log(__a)
#define MY_LCC
//---------------------------------------------------------------------------
#ifndef MY_LCC
/*--------------------------------------------------------------------------*/
double __declspec(naked) myln_sf(double __a){// asm for VC & BC
// 10M Rand data using 1.36 s, log(a) is 1.61 s
asm{
fld dM_LN2 //using fldln2, it is 1.49 s !!!
fld
fyl2x
ret
};
}
#else
/*----------------------------------------------------------------------------------------------*/
double __declspec(naked) myln_sf(double __a){// asm for LCC
_asm("fldl %dM_LN2");
//_asm("fldln2"); // log(2), for speed, intel recommand using a real LN2 number instead of fldln2 which load internal LN2
_asm("fldl 4(%esp)"); // a : log(2)
_asm("fyl2x"); // log(a)
_asm("ret");
}
#endif
/*----------------------------------------------------------------------------------------------*/
#ifndef MY_LCC
void _exportmexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
#else
voidmexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
#endif
{
double *pdAddr;
int i,j;
const int *pnAddr;
if (1!=nrhs) mexErrMsgTxt("Using: myLog( aInOut );Here, aInOut is both input and out put array, for speed up and less memory;");
if(mxDOUBLE_CLASS != mxGetClassID(prhs)) mexErrMsgTxt("Using: myLog( aInOut );aInOut should be a double array;");
pnAddr= mxGetDimensions(prhs);
for(j=mxGetNumberOfDimensions(prhs)-(i=1); j>=0;j--)i*=pnAddr;
for(pdAddr = (double *)mxGetData (prhs);i>0;i--, pdAddr++)*pdAddr = myln_sf(*pdAddr);
}
/*----------------------------------------------------------------------------------------------*/
注，以上默认的是还是LCC，如果是用VC 或者BC编译时请remark "#define MY_LCC"
因为masm和lcc的语法不同，这里给出了两个版本的汇编。
另外，打"//"的语句可以用于VC和BC。dM_LOG2E，dONE 虽然这里没有用到，但是在计算
exp sin cos时就有用了。因为用浮点处理器内部的这些数据速度会比用立即数慢。 写的不错，长知识了 强的～～～～
进来拜一拜
页:
[1]