Matlab 和 C++ （2）

recbio

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 [esp+4]
  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 _export  mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
#else
void  mexFunction(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[0])) mexErrMsgTxt("Using: myLog( aInOut );  aInOut should be a double array;");
pnAddr  = mxGetDimensions(prhs[0]);
for(j=mxGetNumberOfDimensions(prhs[0])-(i=1); j>=0;j--)i*=pnAddr[j];
for(pdAddr = (double *)mxGetData (prhs[0]);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时就有用了。因为用浮点处理器内部的这些数据速度会比用立即数慢。

eisenstange

写的不错，长知识了

rexler_yan

强的～～～～

进来拜一拜