Linux 环境下Fortran程序连接使用Intel MKLPardiso解对称稀疏矩阵

Linux 环境下Fortran程序连接使用Intel MKLPardiso解对称稀疏矩阵

pardiso求解线性方程组[A]{x}={b},其中[A]是对称稀疏阵

<1>Pardiso的Fortran接口：

call pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, perm,nrhs, iparm, msglvl, b, x, error)

简单的参数说明：

pt：长度为64的数组

INTEGER*4 à32位系统

INTEGER*8 à64位系统

说明：内部指针，调用pardiso之前必须置0，之后不要更改，否则会出现内存溢出错误。

maxfct：INTEGER

大多数情况设1

mnum：INTEGER

大多数情况设1

mtype：INTEGER

说明：矩阵形式

结构有限元常用：2 实对称正定矩阵

-2 实对称不定矩阵

流体计算可能会用到 11实数不对称矩阵

phase：INTEGER

一般设为13，analysis->numericalfactorization->solve->iterativerefinement

n: INTEGER

方程阶数

a: REAL/COMPLEX

CSR(compressed sparse row)格式存储的系数矩阵

P.S. CSR格式简介：

设有如下对称矩阵：

CSR格式只存储矩阵上三角中的非0元素，则

a（i）={1，-1，-3，5，4，6，4，7，-5}   

索引ia（i）表示第i行对角线元素在数组a中的位置及对角线后面一个位置，则

ia（i）={1，4，5，8，9，10} ,最后一项ia(n+1)为非零元素个数+1，由于采用上三角存储，非零元素指上三角中的非零元素。

索引ja（i）表示上三角每个非0元素的列标，则

ja（i）={1，2，4，2，3，4，5，4，5}

ia:INTEGER 数组

说明：如上

ja:INTEGER 数组

说明：如上

perm: INTEGER

一般不用，给个虚参

nrhs: INTEGER

方程右边b的列数

iparm: INTEGER 长度64的数组

简单的使用，设iparm（1）=0，则iparm（2），iparm（4）到iparm（64）都将置为默认。iparm（3）为并行计算线程数，一般设为MKL_NUM_THREADS

maglvl: INTEGER

设为1则在屏幕上显示计算统计信息

b: REAL/COMPLEX 数组

方程右边的向量b

x: REAL/COMPLEX 数组

解出的向量x存于这个数组

error:INTEGER

错误代码，有错则输出一个代码，以便查错。

关于pardiso更详细的使用，请参照Intel MKL使用手册，或者登陆http://www.pardiso-project.org

<2>pardiso的和自编程序的链接

准备步骤：确保Linux已经安装Intel Fortran 编译器，和Intel MKL函数库。

本例环境：RedHat Linux Enterprise 5，Fortran编译器及MKL版本10.0

用贯了windows系统，没有图形怎么行，下载vnc_viewer在windows客户端上连接LinuxServer！（过程省略。。。）

正题：

第一步：linux下打开文本编辑器，写一个名为pardiso_sym.f的小程序（本例使用f77固定格式），

如图

具体程序和简单说明如下：

!..Modify from the given example from www.pardiso-project.org

PROGRAM main

IMPLICIT NONE

! 变量定义

INTEGER*8 ipt(64)

INTEGER i,idum

INTEGER maxfct, mnum, mtype, phase, n, nrhs, error, msglvl

INTEGER iparm(64)

Integer mkl_get_max_threads

external mkl_get_max_threads

REAL*8 ddum

integer,allocatable::ia(:),ja(:)

real*8,allocatable::a(:),b(:),x(:)

!为动态数组分配空间,尝试大数组9千万个元素，尝试过系统最大可分配2G

!内存给一个动态数组

!如果数组设为静态，则可以编译通过，但运行时显示Segmentation fault。

WRITE(*,*) 'allocate start'

allocate(ia(90000000))

allocate(ja(90000000))

allocate(a(90000000))

allocate(b(90000000))

allocate(x(90000000))

WRITE(*,*) 'allocate finish'

!各数组赋初值

do i=1,90000000

ia(i)=i

enddo

ia(1:9)=(/1,5,8,10,12,15,17,18,19/)

!注意：如果使用Fortran隐含规则赋初值，如：ia=(/1,5,8,10,12,15,17,18,19,(i,i=10,90000000)/)

!运行时出现Segmentationfault，原因可能是隐含规则赋值调用系统堆栈，系统堆栈默认大小是8192KB，数组长度小则无此问题。

WRITE(*,*) 'ia finish'

do i=1,90000000

ja(i)=i

enddo

ja(1:18)=(/1,3,6,7,2,3,5,3,8,4,7,5,6,7,6,8,7,8/)

WRITE(*,*) 'ja finish'

do i=1,90000000

a(i)=i

enddo

a(1:18)=(/7.d0,1.d0,2.d0,7.d0,-4.d0,8.d0,2.d0,1.d0,5.d0,

1 7.d0,9.d0,5.d0,1.d0,5.d0,-1.d0,5.d0,11.d0,5.d0/)

WRITE(*,*) 'a finish'

n=8 !方程数，虽然数组开的较大，但我们只解8阶方程

nrhs=1

maxfct=1

mnum=1

do i = 1, 64

iparm(i) = 0

end do

iparm(1) = 0 ! iparm使用默认值

iparm(3) = mkl_get_max_threads()

! iparm（3）并行线程数由mkl获得，一般有几个cpu就有几个线程

error = 0 ! 错误信息初值置0

msglvl = 1 ! 显示统计信息

mtype = -2 !对称不定矩阵

do i = 1, 64

ipt(i) = 0

end do

phase = 13 !做全套分析

do i = 1, n

b(i) = 1.d0

end do

WRITE(*,*) 'init finish'

CALL pardiso (ipt, maxfct, mnum, mtype, phase, n, a, ia, ja,

1 idum, nrhs, iparm, msglvl, b, x, error)

WRITE(*,*) 'Solve completed ... '

WRITE(*,*) 'The solution of the system is '

DO i = 1, n

WRITE(*,*) ' x(',i,') = ', x(i)

END DO

WRITE(*,*) 'after call pardiso,a(1:18)=',a(1:18)

C.. 结束及释放内存

phase = -1 ! release internal memory

CALL pardiso (ipt, maxfct, mnum, mtype, phase, n, ddum, idum,

1 idum,idum, nrhs, iparm, msglvl, ddum, ddum, error)

!不要忘记释放动态数组，虽然聪明的系统在你忘了的时候会帮你完成这个任务

deallocate(ia)

deallocate(ja)

deallocate(a)

deallocate(b)

deallocate(x)

END

第二步：打开文本编辑器写一个makefile文件，如图

这里只介绍一下链接pardiso的那一段，就是“ifort –o“的那一段。

Intel Fortran的编译和链接同样由ifort指令完成，

ifort –o pardiso${OBJECTS}表示链接目标文件${OBJECTS}并输出一个名为“pardiso“的可执行程序，紧接着后面-L表示链接库的位置，-I表示include文件或头文件的位置，如图

我的安装的库文件夹位于/ISUM/ProgramFiles/mkl/10.0.3.020/lib/32,（lib文件夹里有对应于32位系统和64位系统的文件夹，这里使用32位编译）同理写出头文件夹的绝对路径。

打开lib/32文件夹，发现里面有很多库文件，以.a结尾的是静态库，.so结尾的是动态库，那么多库，哪些库是pardiso需要的？我也是第一次，还是看看mkluserguide怎么说。

她在mkl安装目录的DOC文件夹里。

userguide给出了我们使用MKL库之前所需明确的东西,我把自己的例子与之对应

①目标机器的架构-->64位, 实际使用了32位编译链接

②要解决的数学问题-->稀疏求解问题,选用pardiso

③编程语言-->fortran

④线程模型-->多线程并行计算,使用intel的openmprun-time库(libguide)

⑤链接形式-->静态链接

⑥是否使用了MPI-->没有使用

之后根据她给的例子使用layered pure link model选择如下库(32位编译):

libmkl_intel.a 使用intel编译器的接口库

libmkl_intel_thread.a 使用intel编译器的并行驱动库

libmkl_solver.a 包含稀疏求解器的库

libmkl_core.a mkl核心库

libguide.so intel自有openmp库

libpthread.so linux系统库

注意:

libguide.so,libpthread.so是动态链接,因为intel强烈建议,即使其他库是静态链接,这两个库也最好动态链接,那么就按她说的办.

P.S.关于intel 层模型的概念内容较多,这里不介绍

如果使用layered pure link model的话,我们需要把libmkl_intel.a,libmkl_intel_thread.a, libmkl_solver.a,libmkl_core.a用一个组符号把他们包起来,即-Wl--start-group<要包含的库的绝对路径> -Wl--end-group

最后不要忘了动态链接libguide和系统库，动态链接的话，把lib改成-l，即libguide改为-lguide，系统库动态链接-lpthread，系统库请放在最后。到这里我们链接mkl到自己编的小程序的链接makefile就完成了。

第三步:现在打开你的终端进入你的编程目录,终端输入make命令开始编译链接:

执行完后,终端显示编译成功,文件夹下出现了,一个名为pardiso的可执行文件.

终端中输入命令./pardiso得到结果

 

下面看个简单例子:

解一个8960阶稀疏矩阵的例子,8个cpu

用skyline方法求解(只用到一个cpu):17秒左右

pardiso并行解她只需要0.1秒