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}
索引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秒
扫一扫
394
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
