找回密码
 注册
Simdroid-非首页
查看: 105|回复: 0

[其他] Eigen入门教程(7)fft和外部库的调用

[复制链接]
发表于 2013-1-11 19:36:18 | 显示全部楼层 |阅读模式 来自 江苏无锡
本帖最后由 myleader 于 2013-1-11 19:57 编辑

fft是对一个数组进行运算,而不是一个矩阵,所以所谓的对一个矩阵进行fft运算其实是对矩阵的每一行或每一列进行fft运算。Eigen官方并没有提供fft功能,不过由于这个功能很常用,所以有热心网友在2009年制作了fft的封装,可以和Eigen一起使用,并成为unsupported的一部分,不过这个模块只是挂着Eigen的名头其实使用起来和Matrix类没有太多关系。甚至Eigen在2011年更新到3版本时,内部接口出现了大幅调整,而fft接口都可以保持不变。

来看一段样例代码

默认的fft使用方法

  1. #include <unsupported/Eigen/FFT>

  2. ...
  3. Eigen::FFT<float> fft;

  4. std::vector<float> timevec = MakeMyData();
  5. std::vector<std::complex<float> > freqvec;

  6. fft.fwd( freqvec,timevec);
  7. // manipulate freqvec
  8. fft.inv( timevec,freqvec);

  9. // stored "plans" get destroyed with fft destructor
复制代码

首先创建了Eigen::FFT类,然后用这个处理一个std::vector类,可以进行变换或者反变换。

现在大家明白楼主为什么说这个模块只是挂着一个Eigen的名头罢了吧,因为它处理的是std::vector类。而且这项功能也不是Eigen自带的代码,而是调用外部库实现的fft算法。

Eigen默认使用的kissfft软件包内部是这样一种机制:先创建一个plan,然后根据输入数据的长度对plan进行优化,然后用这个plan去处理其他的输入数据。这就意味着第一次进行fft时,由于有了创建plan这个额外开销,所以会比较慢;而如果继续使用这个plan进行后面数据的fft,由于省却了创建plan的开销,所以会比较快。

所以如果每次进行fft都要重新创建一个Eigen::FFT对象,那么原有对象中的plan会被销毁,这样每次进行fft都要创建plan,效率就低了。

如果想要对一个矩阵进行逐列fft,那就要使用Map把Matrix类中的数据逐列拷贝到一个std::vector类中,然后对这个vector进行fft,最后再用Map把结果复制回目标Matrix类中。

  1. int m=all.rows(), n=all.cols();
  2. // fft on each column
  3. MatrixXcd result(m,n);
  4. Eigen::FFT<double> fft;
  5. std::vector<double> timevec(m);
  6. std::vector<std::complex<double> > freqvec;
  7. VectorXd tmp_vec_d;
  8. VectorXcd tmp_vec_complex;
  9. for(int n1=0;n1<n;++n1)
  10. {
  11.     tmp_vec_d=all.col(n1);
  12.     Map<VectorXd> (timevec.data(),m,1)=tmp_vec_d;
  13.     fft.fwd(freqvec,timevec);
  14.     Map<VectorXcd> tmp_vec_complex(freqvec.data(),freqvec.size(),1);
  15.     result.col(n1)=tmp_vec_complex;
  16. }
复制代码

例子中使用的输入是实数,如果输入是复数,也没有问题。Eigen是C++写成的,会根据输入数据的类型自动进行重载。

使用其他更快的fft库

有很多软件包提供fft功能,其中最著名的是fftw,它以速度快而著称,不过fftw是严格的GPL协议的,具有开源传染性,所以如果你希望开发闭源软件,那就必须额外买一份fftw的商业授权。Eigen默认使用的是kissfft,这个软件包是BSD协议的,不具有开源传染性,可以用于开发闭源软件,不过它的速度自然是比不上fftw的。另外Intel在其自己的数学库MKL中也提供了fft功能。

实际上Eigen的fft是通过调用外部库来实现的,默认情况下调用的是kissfft,如果你在代码中不加任何声明,其代码中就会
  1. #include<kiss_fft.h>
复制代码

这是不需要我们干预的,Eigen软件包内部已经处理好了,也可以直接用于开发闭源软件。

可使如果你已经购买了fftw或者MKL的商业授权并且希望使用的话,那可以通过修改C++预处理宏来实现调用,只要在文件的头部加上

  1. #define EIGEN_FFTW_DEFAULT
复制代码

要注意这个宏定义必须要在包含Eigen的头文件之前。

Eigen就会调整其内部的头文件和库文件的定义使其指向fftw的头文件和库文件。然后你只要在编译时把fftw.h和libfftw.a放在编译器能找到的目录,然后在编译时添加参数-lfftw就可以了。当然你也可以加上头文件卫士以避免可能的冲突。

kissfft本身仅包含两个文件kiss_fft.h和kiss_fft.c,不需要库文件,所以在编译是不必指定额外的链接库。

使用Intel MKL


通过fft的例子,可以发现在我们使用Eigen调用外部库时,只要代码里定义一些宏就可以了,编译是加上链接外部库的参数就可以了,这大大简化了软件开发的工作,因为我们只要关注前端Eigen提供的接口就可以了,而不必在更换后端计算库时重新编写大量的代码。

在线性代数领域Intel MKL是公认的速度非常快的库。Eigen提供了对Intel MKL的封装,使用起来也很简单,只要定义宏就可以了。更妙的是,Eigen在某些功能上比Intel MKL更好,所以你可以只调用一部分MKL的功能,而其他部分保留Eigen的后端。这些宏包括:

EIGEN_USE_BLAS 使用Intel MKL中BLAS level2和level3的功能
EIGEN_USE_LAPACKE 调用Intel MKL中LAPACKE的功能
EIGEN_USE_LAPACKE_STRICT 和 EIGEN_USE_LAPACKE类似,不过一些鲁棒性较差的函数不会被调用。最主要的问题在于JacobiSVD就是教程(4)里介绍的最小二乘法矩阵分解。
EIGEN_USE_MKL_VML 调用Intel MKL中的VML功能。
EIGEN_USE_MKL_ALL 同时包含EIGEN_USE_BLAS、EIGEN_USE_LAPACKE和EIGEN_USE_MKL_VML

针对不同的应用,Eigen会调用MKL中不同的函数。对于不同的应用所调用的函数,这里简单介绍一下:
运算
EIGEN_USE_BLAS 教程(3)中介绍的矩阵算数运算,如转置、伴随矩阵、三角阵、加减乘除等
EIGEN_USE_LAPACKE 教程(4)中介绍的矩阵分解,如LU、QR、LLT、LLDT、SVD等
EIGEN_USE_MKL_VML 教程(3)中介绍的逐元素运算,如sin、cos、tan、cot、asin、acos、atan、exp、ln、pow等

可以根据需要调用必须的部分。实际上如果你非要用AVX来完成矩阵的算术运算理论上也是可以的,不过Eigen内部类和函数的封装中使用了MKL中的SSE指令,而不是AVX指令,当然Eigen自己的代码中使用的也是SSE。所以其实MKL不会更快。

MKL后端和Eigen后端的对比

MKL是Intel搞出来的,用在Intel的CPU上当然速度杠杠的,不过如果你使用的是非Intel的CPU,就要多考虑一下了。Eigen目前内部的矢量化代码支持Arm的NEON指令集,还支持IBM PowerPC的AltiVec指令集,所以如果你是用这两种CPU的话,仍然能得到Eigen的加速。不过目前还不支持MIPS的矢量指令集

另外MKL中AVX的部分目前仍然只能单线程,这也是一个硬伤,如果你的问题更适合多线程的话,不如使用SSE多线程计算速度还要快一点。

所以说不能迷信MKL,要根据你自己的实际题目来选择后端。

EIGEN_USE_BLAS调用了BLAS level2和level3的功能,其实这个东西意义不大,因为Eigen的速度已经很快了,如果调用MKL还要多一层动态库的载入,数据量小时反倒不划算。

EIGEN_USE_LAPACKE,其实LAPAKE主要包含的就是教程(4)里面提供的几种分解计算,而Eigen的分解计算也已经非常快了,不过由于都是大数据处理,所以动态库载入的消耗可以略去不计。如果你确认你的算法中不会用到JacobiSVD分解,用EIGEN_USE_LAPACKE_STRICT也未尝不可,自己心里放心点。

EIGEN_USE_MKL_VML调用Intel MKL中的VML功能。Intel从Sandy Bridge推出AVX指令集,AVX在矢量计算中表现出了极强的力量,速度比SSE快了好多。如果你的CPU支持AVX指令集,而你的题目也适用AVX的话,那么楼主推荐。

不过还有一个问题就是截止2012年底MKL中的AVX指令集功能仍然不支持多线程。据OpenBLAS的测试,如果只能单线程的话带AVX的OpenBLAS的速度甚至比MKL还要快。也就是说Intel在这个指令集上挖掘出来的潜能还没有挖光,但是从另外一个角度看来,如果你的CPU数量巨大,足以弥补不能使用AVX指令集带来的遗憾的话,从整体的计算速度考虑,还是退回到SSE多线程方案比较好。

综合起来,如果你的数据量比较小,那么楼主不推荐使用Intel MKL作后端,如果你的数据量中等,单线程AVX指令集最合适的话(主要是三角函数、乘方和对数的运算),而且你的硬件还支持AVX指令集,那么楼主推荐。不管是巨量数据(可以用到很多线程)还是说你的CPU不支持AVX,楼主都不推荐MKL后端。

虽然OpenBLAS也不差,不过Eigen目前还不能支持OpenBLAS后端。

其他外部库

其实楼主在教程(5)中就已经介绍了一些外部库的使用,主要是稀疏矩阵,对于那些库的调用,其实不需要额外的设置,因为当你创建矩阵分解器时就已经在include相关的头文件了,如果编译器找不到对应的头文件,编译就会出错。所以那些模块的使用只要让编译器能找到对应的头文件和库文件就好了。至于那些软件包的安装,楼主如果有空再给大家介绍。

其实官方正式模块在教程(6)中就已经介绍完了,接下来的非正式模块中,比较常用的除了fft以外也就只有多项式Polynominal模块用得比较多,不过非正式模块是没有教程的,只能看代码学习,楼主没怎么看懂,所以恐怕一时半会不会写了。暂时就先到这里,如果有任何进展,楼主觉得值得一写,那就再继续,如果另有高人,也欢迎来介绍。

您需要登录后才可以回帖 登录 | 注册

本版积分规则

Archiver|小黑屋|联系我们|仿真互动网 ( 京ICP备15048925号-7 )

GMT+8, 2024-3-29 08:12 , Processed in 0.027298 second(s), 11 queries , Gzip On, MemCache On.

Powered by Discuz! X3.5 Licensed

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表