快速多极子方法并行的技术

s( x j ) i ( x j xi ),{s j } [ ji ]{ j }
i 1
N
• • • •
以上求和要求O(MN)运算复杂度 对给定的精度，FMM可以获得O(M+N)运算复杂度 可以加速matix-vector products ，使O(N2)变为O(NlogN) 加速线性系统求解，如果用迭代方法，k步收敛，每步用矩阵矢量 相乘，使计算复杂度由O(N3)变为O(kNlogN)
i
N
li ρ (r ) / 2 A , ρ (r) r ri , if r Ti ; fi (r ) (r ) / 2 Ai , ρi (r) ri r, if r Ti . li ρi N
+
i 1
•
z j v , m 1, 2, , N . MOM 离散 z f (r) L(f ) (1 )[f (r) / 2 n K (f )] dS ,

每个盒子在l=0,..,L层的指标n=Number=0,1,…,2ld-1.
E1 (n, l ) R d 表示盒子(n, l )的空间区域; E2 (n, l ) R d 表示盒子(n, l )与其相邻{( Neighbor (n, l ), l )}的空间区域; E3 (n, l ) E1 (0,0) \ E2 (n, l )表示盒子(n, l )及其相邻{( Neighbor (n, l ), l )}的外部空间区域； E4 (n, l ) E2 ( Parent (n), l 1) \ E2 (n, l )表示父盒子( Parent (n), l 1)及其相邻{Neighbor(Parent (n), l 1), l 1}的空间区域；
3
FMM in Computational Electromagnetics
• EFIE
t L(J ) t Ei (r) , r S , L(J) jk I / k
S 2
g (r,r') J(r')dS ';
• MFIE
t J (r ) / 2 t n K (J ) t n Hi (r) , K (J ) J (r') g (r,r')dS ' ;
国家973项目——高性能科学计算研究 大规模并行计算研究
快速多极子方法的并行技术
冯仰德 王武 迟学斌
中科院计算机网络信息中心 超级计算中心 ydfeng@sccas.cn 2007年2月5日
1
纲要
• FMM • Data Structures • Parallelization
2
纲要
• FMM
• Data Structures • Parallelization
8
FMM:Application
• Molecular and Stellar dynamics —computation of force fields and dynamics • Solution of acoustical scattering problems —Helmholtz Equation • Electromagnetic Wave Scattering — Maxwell’s Equations • Fluid Mechanics:Potential flow,vertex flow — Laplace/Poisson Equations
xi E1 ( n ,l )
ui (y, x i ), (2n ,l ) (y ) ui (y, xi ), (4n ,l ) (y )
xi E2 ( n ,l )
ui ( y, x i ), ui ( y, x i )
xi E3 ( n ,l )

xi E4 ( n ,l )
注意到l<z时，函数jl(z)和hl(2)(z) 幅值大致保持为常数； l>z时，函数jl(z)衰减非常快， 而hl(2)(z)递增非常快。当d<<r时，上式在保证精度的情况下截断。则上式可以写为：
L e jk |r+d| jk (1)l (2l 1) jl (kd )hl(2) (kr ) Pl (d r) |r +d| l 1
Based On: ◆ Hierarchical Grouping ◆ Data Structure
18
Multilevel Multipole Operators
Finest - 1 Level
M2L
M2M
M2M
M2L
L2L
L2L
Finest Level
M2L
Up Tree
Across Tree
7
Fast Multipole Methods(FMM)
• Introduced by Rokhlin &Greengard in 1987 • Called one of the 10 most significant advances in computing of 20th century • Speeds up matix-vector products (sums) of a particular type
M2M , M2L , L2L： 聚集 – 转移 – 发散 M: 多极子展开 L: 局部展开
13
FMM
(n) 1 (y )
xi E1 ( n )

ui ( y , x i ), ui ( y , x i ), ui ( y , x i )
(2n ) ( y )
(n) 3 (y )
M2L
16
Step3
L2L
Summation
17
MLFMM Algorithm
• Upward Pass: merge scattering matrices • Downward Pass: construct splitting and exchange matrices • Final Summation
Pulse
• Basis Functions Interact through the Green’s Function
f (r )
G r, r
t (r )
• Generates a Dense Method of Moments Matrix
ZnnI n1 Vn1
Z j, i f j , G r, r ' fi r ' ds
Down Tree
19
MLFMM Algorithm —Upward Pass
• Step1: 在最细的层盒子求解远场展开系数，xi∈E1(n,l), 1(n,l ) (y) 得到C(n,l)或 ，这也可以用于 xi∈E3(n,l) • Step2:对于l=L-1,…,2,从step1得到值进行递归得到。 同样适合xi∈E3(n,l) 结果：得到分层组聚集系数
xi E2 ( n )

由于E2(n)和E3(n)是互补的，因此对任意的n,
xi E3 ( n )

(y)
xi E2 ( n )E3 ( n )

( n) ui (y, xi ) (2n) (y) 3 ( y)
14
FMM Algoritห้องสมุดไป่ตู้m
Step1 M2M
15
Step2
NeighborsAll(n,l), Neighbors(X;n,l)
23
Grouping

由于E2(n,l)和E3(n,l) 是互补的，因此对 任意的n,l
n,l ) ( n,l ) (y) ( (y) 3 (y ) 2
( n ,l ) 1 (y ) ( n ,l ) 3 (y )
i
Dip e jk k(r rm ) [ (I - k p k p ) fi (r) (1 )k p fi (r) n]dS
Tp TL (krmm ' , k p r mm ' ), L kd c ln(kd ), K 2 L2 ; A pj e jk k (r 'rm ' ) (I - k p k p ) f j (r ')dS ' ;
S
6
线性系统：
Mx=s M是N×N矩阵，x、s是N矢量 ●Direct solution (Gauss elimination,LU Decomposition,SVD,…) 空间复杂度为O(N2) ，需要O(N3)次 运算； ●Iterative methods,空间复杂度仍为O(N2)，如果K(k<<N)步收 敛，每步需要的矩阵乘向量的运算为O(N2)； — Moore’Law:processor speed doubles every 18 months — a million variable,16 generations of Moore’Law before a O(N2) algorithm was comparable with a O(N) algorithm — 1 GB RAM =10243=1,073,741,824 bytes => largest N =32,768 ●Finding:快速矩阵乘向量的算法(NlogN)； 并行实施。
mn n m n=1
mn S m n n n i i vm fm (r) E (r) (1 ) n H (r) dS. S
5
Method of Moments
• Surface is Discretized into Patches (Basis Functions)
S
• CFIE CFIE • Green函数
EFIE (1 ) MFIE
jk |r -r '|
g (r, r ') e
/ 4 | r - r ' |
4
积分方程的离散
Rao-Wilson-Glisson
• RWG矢量基函数
i i
J (r ) ji fi (r )
20
MLFMM Algorithm
—Downward Pass
• Step1:l=2,..L递归进行E4
任意一个非空组自身加上其邻接组最多可 有3d个，其父组及父组的邻接组最多可形 成3d×2d，因此次相邻数目最多为 p=3d(2d-1)；三维是189，二维是27，一 维是3。
• Step2:
局部到局部的转移，E3(n,l)和E4(n,l+1) 产生E3(n,l+1)
9
FMM: Fundament
• 格林函数的加法定理
e jk (1)l (2l 1) jl (kd )hl(2) (kr ) Pl (d r) |r +d| l 1
jk |r+d|
jl_—第一类球面Bessel函数 hl(2) —第二类球面Hankel函数 Pl —Legendre多项式
j
12
FMM
{ y} { y}Near [ D][T ][ A]{x}
{x1} [ A]{x}
—— 聚集过程，将基函数聚集成平面波函数，其结果表示K个平面波
{x2 } [T ]{x1} ——转移过程，将得到的{x1} 平移到另一个盒子的中心，其结果表
示该盒子中心的K个平面波
{ y1} [ D]{x2} ——发散过程，将得到的{x2}发散到盒子中的基函数。
L=kd+cln(kd+ )
Kd-源点到观察点的最大半径 c-是一个依赖希望精度的常数 =1 最小的相对误差小于0.1 =5 相对误差小于10-6 =10 准确到双精度
• jlpl平面波展开
4 ( j ) n jn (kd ) P n (d r )

Sa
2 e jk k d P ( k r ) d k n
10
Fast Multipole Basics
直接求解： 分解：
{cm}
复杂度：O(MN)
复杂度：O(pN+pM)
11
电磁场积分方程的多极子展开
• FMM形式的矩阵向量乘积
z
n 1
N
ij
jj
m 'Gn ( m ) jGm '
z
近场部分
ij
jj
远场部分
k 2 K w D [ Tp ( A pj j j )] 2 p ip 16 p 1 m 'Gn ( m ) jGm '
结果：可以得到各分层组的转移系数
21
Key Words
• • • • •
空间多层组划分 Morton编号 相邻组的作用 远场组的上聚 次相邻组中心 的转移 • 远场组的下推
22
Grouping
● 每个盒子在第l层的指标号数为n,那么任意盒子的指标为(n,l); ● 需要构建的函数，如Parent(n), ChildrenAll(n), Children(X;n,l),