24小时热门版块排行榜

返回列表

本帖产生 1 个程序强帖，点击这里进行查看

coolrainbow

木虫 (著名写手)

未来国家冻凉

程序强帖: 2
应助: 1 (幼儿园)
贵宾: 0.1
金币: 3979.5
散金: 10
红花: 50
帖子: 1137
在线: 272.3小时
虫号: 239560
注册: 2006-04-07
性别: GG
专业: 理论和计算化学

[交流] 【分享】MKL程序编译与连接：Lapack篇已有5人参与

原帖转自我的空间：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[j]的形式存储于计算机中。这是最简单的方法，但是对稀疏矩阵或元素有一定相关性的矩阵效率较低

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，这个算矩阵乘法的程序代码：

CODE:

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：

CODE:

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，汗颜....

[ Last edited by coolrainbow on 2010-11-18 at 15:21 ]

回复此楼