MKL程序编译与连接：Lapack篇

 原帖转自我的空间：http://hi.baidu.com/coolrainbow/blog/index/1

经验表明，第一次做MKL程序编译时，大多数人都会走很多弯路，编译几个小时也不成功。其实这个原因很简单，就是：懒于阅读技术文档！！！！！！事实上，本人当时就是这样，看到3000多页的MKL手册头都大了，于是在编译时就“跟着感觉走”，弄得找不到库或者函数错误。其实，花一个小时读下文档，绝对比自己瞎折腾要强的多，这里把编译MKL的一些经验与大家分享下，作为快速入门。真正深入的话，还是那句：读文档！

1 MKL的环境变量

安装好MKL后，需要设置一些环境变量，这样才能找到所需要的库，这可以通过/opt/intel /Compiler/11.1/064/mkl/lib/tools/enviroments/mklvars{your-architecutre}. {sh|csh}实现。如果需要的话，可以加入到/etc/profile或你的.bashrc中。在程序编译或运行时，如果发生can not find libXXX之类，记得导出相应的LD_LIBRARY_PATH。如：

can not find libmkl_intel_thread: cannot open shared object file‎...

如果你的这个库位于/opt/intel/Compiler/11.1/064/mkl/lib/em64t，那么：

export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/opt/intel/Compiler/11.1/064/mkl/lib/em64t

即可。

2 函数选择：Lapack

MKL提供了3000多页的技术文档，对每一个函数进行了详细的介绍。这里仅对Lapack系列做一些介绍。

初学Lapack的人可能会问：我怎么知道哪个函数是做向量点积的，哪个函数是做奇异值分解的。问的最多的是：有没有什么书可以参考的？答案：至少我没见过哪个书像讲MFC或PHP那样详细的把Lapack的每一个函数讲解的书。那么唯一的途径就是阅读技术文档。不要头疼，那个文档有英文版和日语版（shit）的，没有中文版的。好好学英文吧。

Lapack中的函数的名称都是XYYZZZ形式（BLAS类似）的，

X     精度。 如
s     单精度（对应于Fortran中的real(kink = 4)，C中的float）
d     双精度（对应于Fortran中的real(kink = 8)，C中的double）
c     单精度复数   z      双精度复数

YY   对应的矩阵类型和存储方式

ge    一般的矩阵，存储方式full
所谓full就是说，矩阵A的元素aij以数组a(i, j)或a的形式存储于计算机中。这是最简单的方法，但是对稀疏矩阵或元素有一定相关性的矩阵效率较低

sy     对称矩阵，以full方式存储
sp     对称矩阵，以pack方式存储
pack方式对symmetric matrix采用了一维数组优化存储的方法，可以减小内存的使用。有两种存储方法，在Lapack程序中分别用'U' 'L'表示，对于'U'，简单的说：
aij --->    a(i+j*(j-1/2)) where i < j
'V‘类似。

还有其他的类型，请参考技术文档。

ZZZ 任务类型

就是你要干什么。如dot表示点积，trd表示三对角化等等

知道了这个命名规则后，找函数就方便多了。

3 函数参数：Lapack

Lapack中的每一个函数都有很长的参数，常常把初学者吓倒。并不是说对矩阵A，B求乘法，只要简单的输入A，B就可以了。还有一些辅助空间，任务类型之类的，都要输入，如矩阵乘法：

call dgemm('n','n',N,M,K,a,x,N,y,K,b,z,N)

执行矩阵乘法z=x*y，但还有一些character*1 的'n','n'以及integer的N,K等等,这些含义一定要仔细弄清，不然可能会使任务执行错误！

4 编译连接

这是最头痛的问题。MKL的连接选项比较长，最好使用make工具，我常用的一个模板如下

mkllib=/opt/intel/Compiler/11.1/064/mkl/lib/em64t
mklinc=/opt/intel/Compiler/11.1/064/mkl/include
foo: foo.o
           ifort -o foo foo.o -I$(mklinc) -L$(mkllib) -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -lmkl_lapack95_lp64 -liomp5 -lpthread

在连接MKL程序时，要在每一个“LAYER”中选择一个库，
1 interface layer：提供程序接口，即libmkl_intel{lp64}，其中lp64是对em64t或ia64版本。
2 threading layer：线程接口。一般都是libmkl_intel_thread。注意，有了它后，一定要加上iomp5和pthread库，不然会出现连接错误。
3 computational layer，计算接口，即libmkl_core和你所需要的，如lapack的libmkl_lapack95_lp64或者线性方程组的libmkl_solver_lp64等等。
4 RTL接口，即pthread等等。

这样才可以连接成功。

实例：

现在以一个实例完成。想算一个矩阵乘法，矩阵用双精度，一般存储模式，则选用dge???函数，查文档得知?gemm做矩阵乘法，于是选用dgemm，这个算矩阵乘法的程序代码：


program main
    integer,parameter::N=3000
    integer,parameter::M=2000
    integer,parameter::K=4000
    real(kind=8) a,b
    real(kind=8),pointer::x(:,:)
    real(kind=8),pointer::y(:,:)
    real(kind=8),pointer::z(:,:)
        a=1.0
        b=0.0
        allocate(x(N,K))
        allocate(y(K,M))
        allocate(z(N,M))
        x=10.0
        y=20.0
        z=0.0
        write (*,*) "Performing matrix multiplication..."
        call dgemm('n','n',N,M,K,a,x,N,y,K,b,z,N)
        write (*,*) "DONE!"
   deallocate(x)
   deallocate(y)
   deallocate(z)
end


命名为dge.f90，写Makefile：


mkllib=/opt/intel/Compiler/11.1/064/mkl/lib/em64t
mklinc=/opt/intel/Compiler/11.1/064/mkl/include
FCCFLAG= -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -lmkl_lapack95_lp64 -liomp5 -lpthread
FCC = ifort
dge: dge.o
        $FCC -o dge dge.o -I$(mklinc) -L$(mkllib) $(FCCFLAG) -g
dge.o: dge.f90
        $FCC dge.f90 -c -g


然后执行make，好，一个叫dge的可执行程序出现了，在我的机器上运行时间为1.457s。用我自己随手写的一个矩阵乘法程序一算：11min34s，汗颜....

http://emuch.net/html/201012/2608578.html