to AztecOO Department of Science and Engineering Computing School of Mathematics School Peking University October 9, 2007

Size: px

Start display at page:

Download "to AztecOO Department of Science and Engineering Computing School of Mathematics School Peking University October 9, 2007"

剧姬
5 years ago
Views:

1 to AztecOO Department of Science and Engineering Computing School of Mathematics School Peking University October 9, 2007

2 to AztecOO Department of Science and Engineering Computing School of Mathematics School Peking University October 9, 2007

3 框架 1 2 Method of Undetermined Coecients BVP 3 Second-Order Initial Value Problems

4 框架 Method of Undetermined Coecients BVP 1 2 Method of Undetermined Coecients BVP 3 Second-Order Initial Value Problems

5 两个隐式格式 Method of Undetermined Coecients BVP 一阶导数的四阶隐式差分格式 1 [(u x) i+1 + 4(u x ) i + (u u i+1 u i 1 x ) i 1 ] = + O( x 4 ) 6 2 x 要求解一个三对角线性方程组才能解的一阶导数的差分近似 (u x ) i 和一阶导数的四阶差分格式比较, 显示格式具有更宽的 stencil. (u x ) i = u i+2 + 8u i+1 8u i 1 + u i 2 12 x 二阶导数的四阶近似 (u xx ) i (u xx ) i + (u xx ) i O( x 4 ) = u i+1 2u i + u i 1 x 2 +O( x 4 ) 要求解一个三对角线性方程组才能解的二阶导数的差分近似 (u xx ) i

6 框架 Method of Undetermined Coecients BVP 1 2 Method of Undetermined Coecients BVP 3 Second-Order Initial Value Problems

7 Method of Undetermined Coecients Method of Undetermined Coecients BVP 5/7 Stencil: c s s s s s i 2 i 1 i s s s s s s s i 3 构造隐式格式 b i 2 a i 1 i i + 1 a i + 1 i + 2 b i + 2 LHS c RHS i + 3 (u x ) i 2 + (u x ) i 1 + (u x ) i + (u x ) i+1 + (u x ) i+2 = c u i+3 u i 3 6 x + b u i+2 u i 2 4 x + a u i+1 u i+2 2 x

8 Method of Undetermined Coecients BVP 3/5Stencil i-1 i i+1 i-2 i-1 i i+1 i+2 = 0 得到三对角矩阵又 c = 0, 得到 3=5 Stencil 格式. a = 2 ( 3 + 2); b = 1 (4 1) 4 LHS RHS

9 一些格式 Method of Undetermined Coecients BVP = 1=4; a = 3=2; b = 0, 3=3 Stencil. O( x 4 ) = 1=3; a = 14=9; b = 1=12, 5=3 Stencil. O( x 6 )

10 The Second order derivative Method of Undetermined Coecients BVP (u xx ) i 1 + (u xx ) i + (u xx ) i+1 = a u i+1 2u i + u i 1 + b u i+2 2u i + u i 2 x 2 4 x 2 一个四阶格式 = 1=10; a = 6=5; b = 0 得到 classical 3=3 Pade scheme. 一个六阶格式 = 2=11; a = 12=11; b = 3=11 得到一个六阶格式

11 框架 Method of Undetermined Coecients BVP 1 2 Method of Undetermined Coecients BVP 3 Second-Order Initial Value Problems

12 Method of Undetermined Coecients BVP Poisson equation u 00 (x) = f (x); x 2 [0; 1]; x(0) = x(1) = 0 二阶显示离散 u i 1 + 2u i u i 1 = x 2 f i ; i = 1; ; N 1 四阶 u i 1 +2u i u i+1 = x2 12 (f i f i + f i+1 ) ; i = 1; ; N 1 作业 : 选取不同的 N, 比较两种格式的精度, 计算数值收敛的阶

13 框架 Second-Order Initial Value Problems 1 2 Method of Undetermined Coecients BVP 3 Second-Order Initial Value Problems

14 问题 : 1st order ODEs Second-Order Initial Value Problems 方程 Euler-Backward Scheme d dt U = F(u; t); U (t = 0) = U 0 U n+1 U n t Crank-Nicolson Scheme U n+1 U n t = F(U n+1 ; t n+1 ) = F(U n+1 ; t n+1 ) + F(U n ; t n ) 2

15 Second-Order Initial Value Problems Third-order Backwards Dierentialtion Method 11U n+1 18U n 9U n 1 + 2U n 2 = 6F n+1 t stiy stable method ( 刚性稳定 ) Third-order Adams-Moulton method U n+1 t U n = 5 12 F n F n 1 12 F n 1

16 框架 Second-Order Initial Value Problems 1 2 Method of Undetermined Coecients BVP 3 Second-Order Initial Value Problems

17 Second-Order Initial Value Problems Second-Order Initial Value Problems 二阶初值问题 U tt (t) = F(t; U (t); U t (t)); U (0) = U 0 ; U t (0) = Z 0 化成一阶问题 U t = Z; Z t = F(t; U (t); Z)

18 二阶问题的 Newmark 方法 Second-Order Initial Value Problems 二阶精度无条件稳定的一个 1 U n+1 = U n + tz n + t 2 F n F n 4 4 Z n+1 = Z n + t 1 F n F n 2 2 二者合并可以写成 U n+2 2U n+1 + U n = t F n F n F n

19 框架 1 2 Method of Undetermined Coecients BVP 3 Second-Order Initial Value Problems

20 怎样选择时间积分方法 I How to March in Time The Rung-Kutta method is the easiest to use and is often a good choice, although somewhat inecient. Low storage and TVB (Total Variation Bounded) schemes are useful in practice. Boundary conditions at intermediate stages may be an issue as they may aect accuracy. Multi-step methods (e.g., Adams-Bashforth) are often used in discretizing convective contributions. They require a lot of memory, but they are ecient. The high-order versions exhibit very low dispersion errors. In an unstable multi-step time integration the error has oscillatory form. Leap-frog schemes are good for problems with complex eigenvalues.

21 怎样选择时间积分方法 II Predictor-corrector methods are good except when it is necessary to change the step size. It is also easy to estimate time errors. They are less stable than implicit. Sti equations require implicit methods. For eciency, stiy stable methods (i.e., backwards dierentiation) are recommended. The truncation error of a k-order stiy stable method is (1=k) compared to the corresponding Adams-Moulton method. Typical instabilities appear in a spontaneous (explosive) way it only takes a few time steps for the solution to blow up!

22 框架 1 2 Method of Undetermined Coecients BVP 3 Second-Order Initial Value Problems

23 框架 1 2 Method of Undetermined Coecients BVP 3 Second-Order Initial Value Problems

24 Serial Thomas 三对角矩阵线性方程组 Ax = q A 是三对角矩阵,LU 分解如下 0 a1 c1 b2 a2 c2 0 b3 a3 c b N a N 1 C A = LU where L = 0 1 l2 1 0 l l N 1 1 C A ; U = 0 d1 u1 d2 u2 0 d3 u d N 1 C A

25 求解三步骤 1 LU 分解 ( 假定多存一些数据, 令 b 1 = 0 and c N = 0 ) l i = b i =d i 1 ; d i = a i l i u i 2 Forward substitution 求解 Ly = q 1 ; u i = c i y 1 = q 1 ; y i = q i l i y i 1 ; i = 2; ; N 3 Backward substitution 求解 Ux = y x N = x N =d N ; x i 1 = y i 1 u i 1 y i ; i = N 1; ; 1

26 Thomas 的递推公式对角元 d i l i = b i =d i 1 ; d i = a i l i u i 1 ; u i = c i ) d i = a i b i u a i 1 id i 1 b i c i 1 = d i 1 d i 利用, 我们有 1 0 d i = 0 1 T i 1 1 where T i = R i R i 1 R 1 T i 1 1 R i = ai b i c i 1 1 0

27 框架 1 2 Method of Undetermined Coecients BVP 3 Second-Order Initial Value Problems

28 transformation I 下面的变换被称作变换 t = A! + B C! + D 在并行计算的课上, 我们要用到变换的 2-by-2 矩阵表示考虑两个变换, 一个是上面给出的形式, 另一个是! = az + b cz + d 两个变换的复合是简化上面的公式, 得到 t = A az+b cz+d + B C az+b cz+d + D t = (Aa + Bc)z + (Ab + Bd) (Ca + Dc)z + (Cb + Dd)

29 transformation II 上面出现的系数刚好是下面两个矩阵的乘积 A B a b = Aa + Bc Ab + Bd C D c d Ca + Dc Cb + Dd 复合运算可以用矩阵乘法来表示这是两个复合, 两个矩阵相乘, 自然可以推广到任意多个变换的系数如果是 A, B, C, D ( 第一个公式 ) 则这个变换可以写成 1 T A B 0 t = 0 1 C T A C D B D! 1! 1

30 transformation III 单位变换是 t =! 可以写成 1 0 t = 0 1 T A C T 1 0 B D 0 1! 1! 1 一个 2-by-2 矩阵实际上和一个变换有着一一对应的关系因此就可以用 2-by-2 的矩阵来表示变换

31 框架 1 2 Method of Undetermined Coecients BVP 3 Second-Order Initial Value Problems

32 前缀扫描 P i ; i = 1; 2; ; N 处理器上分别存放 a i, 要求的计算结果是 b i = ix j=1 例如 : f1; 2; 3; 4; 5; 6; 7; 8g 则计算结果是 f1; 3; 6; 10; 15; 21; 28; 36g a i

33 算法 1: 示意图第一步 :P i 的数加到 P i+1 上, i = 1; 且 i + 1 N 第二步 : 将 P i 的数加到 P i+2 上, i = 1; 且 i + 2 N 第三步 : 将 P i 的数加到 P i+2 2 上, i = 1; 且 i + 4 N

34 算法 2: 示意图 I 改错 : 图中 2 i 应该为 2i, 4 i 应该为 4i,

35 算法 2: 示意图 II 第一步 :P 2i 1 的数加到 P 2i 上, i = 1; ; 2i N 第二步 : 将 P 2 2 i 2 的数加到 P 2 2 i 上, 2 2 i 2 0; 2 2 i N 第三步 : 将 P 2 3 i 2 2 的数加到 P 2 3 i 上, 2 3 i 2 2 0; 2 3 i N 第四步 : 将 P 2 2 i 的数加到 P 2 2 i+2 上, i = 1; ; 2 2 i + 2 N 第五步 : 将 P 2i 的数加到 P 2i+1 上, i = 1; ; 2i + 1 N

36 MPI_Scan 实现上面的运算 I int MPI\_Scan( void *sendbuf, void *recvbuf, int count, MPI\_Datatype datatype, MPI\_Op op, MPI\_Comm comm ); Parameters sendbuf [in] 发送数据开始的地址 recvbuf [out] 接收数据开始的地址 count [in] 发送的长度 datatype

37 MPI_Scan 实现上面的运算 II [in] 数据类型, 例如 MPI_INT op [in] 操作,e.g. MPI_SUM comm [in] communicator, e.g. MPI_COMM_WORLD 上面例子可以使用下面的语句实现 int rank = MPI_Comm_rank(MPI_Comm_WORLD,&rank); int a, b; a = rank +1; MPI_Scan(&a,&b,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);

38 自己定义二目运算 MPI_Op I /** * \file mpiscan.cpp 前缀归约操作使用自己创建的二目运算符 * * ===================================================================================== * * Filename: mpiscan.cpp * * Description: * * Version: 1.0 * Created: 2007 年 10 月 30 日 12 时 49 分 37 秒 CST * Revision: none * Compiler: mpicxx -DMPICH_IGNORE_CXX_SEEK -o mpiscan.exe mpiscan.cpp * * Author: Tiao Lu, tlu@math.pku.edu.cn * Company: School of Mathematical Sciences, Peking University * * ===================================================================================== */

39 自己定义二目运算 MPI_Op II #include <stdio.h> #include "mpi.h" /** 定义个二元运算符 * * invec 和 inoutvec 都是长度为 4 的向量, 分别对应 * 2-by-2 的矩阵 A 和 B, * A[0][0]=invec[0]; A[0][1]=invec[1]; * A[1][0]=invec[2]; A[1][1]=invec[3]; * 这个 operator 给出矩阵的乘积 * B = A*B * invec IN inoutvec OUT length IN

40 自己定义二目运算 MPI_Op III data_type IN */ void oper(void * invec, void * inoutvec, int * length, MPI_Datatype * data_type) { int i; double * a = (double *)invec; double * b = (double *)inoutvec; double b0tmp, b1tmp, b2tmp; b0tmp = a[0]*b[0] + a[1]*b[2]; b1tmp = a[0]*b[1] + a[1]*b[3]; b2tmp = a[2]*b[0] + a[3]*b[2]; b[3] = a[2]*b[1] + a[3]*b[3]; b[0] = b0tmp; b[1] = b1tmp; }

41 自己定义二目运算 MPI_Op IV int main(int argc, char * argv[]) { int rank; MPI_Op op; double buffer[4] = {1, 0, 2, 3}; MPI_Init(&argc, &argv); MPI_Comm_rank(MPI_COMM_WORLD, &rank); //MPI_Op_create: 第一个参数是二元运算函数, // 第二参数非零表示满足交换律, // 第三个参数返回 MPI_Op 运算符 MPI_Op_create(&oper, 0, &op); // 前缀归约使用自己定义的 MP_Op op MPI_Scan(buffer, buffer, 4, MPI_DOUBLE,

42 自己定义二目运算 MPI_Op V op, MPI_COMM_WORLD); MPI_Op_free(&op); MPI_Finalize(); printf("rank %d: {%lf, %lf, %lf, %lf}\n", rank, buffer[0], buffer[1], buffer[2], buffer[3]); return 0; } 这个操作稍微改变一下, 可以用来简化三对角线性方程的并行 Thomas 算法的程序

43 框架 1 2 Method of Undetermined Coecients BVP 3 Second-Order Initial Value Problems

44 I 函数的头部说明头文件和说明 /* * ===================================================================================== * * Filename: thomasp.cpp * * Description: 求解三对角方程组的 Thomas 算法的并行实现 * 使用代替原来的手动操作 * 使用前缀扫描函数可以不使用辅助变量 T * 但是二者都要求使用的 processor 的个数是 \f$ 2^n \f$, * 不然就会报错或者计算结果不正确 * * Version: 2.0 * Created: 2007 年 10 月 30 日 11 时 02 分 14 秒 CST * Revision: none * Compiler: mpicxx -DMPICH_IGNORE_CXX_SEEK -o thomasp.exe thomasp.cpp * * Author: Tiao Lu, tlu@math.pku.edu.cn * Company: School of Mathematical Sciences, Peking University, P. R. China * * ===================================================================================== */ #include <iostream> #include <math.h> #include "mpi.h"

45 II 二目运算符计算两个 2-by-2 矩阵的乘积 using namespace std; /** 定义个二元运算符 * * invec 和 inoutvec 都是长度为 4 的向量, 分别对应 * 2-by-2 的矩阵 A 和 B, * A[0][0]=invec[0]; A[0][1]=invec[1]; * A[1][0]=invec[2]; A[1][1]=invec[3]; * 这个 operator 给出矩阵的乘积 * B = B*A * invec IN inoutvec OUT length IN data_type IN */ void oper(void * invec, void * inoutvec, int * length, MPI_Datatype * data_type) { int i; double * a = (double *)invec; double * b = (double *)inoutvec; double b0tmp, b2tmp; b0tmp = b[0]*a[0] + b[1]*a[2]; b[1] = b[0]*a[1]+b[1]*a[3];

46 III b2tmp = b[2]*a[0] + b[3]*a[2]; b[3] = b[2]*a[1] + b[3]*a[3]; b[0] = b0tmp; b[2] = b2tmp; } 并行 Thomas 算法的函数头写成一个单独的函数, 可以重用 \** 并行追赶法. * * 三对角矩阵是 (b, a, c) * mynode IN 进程号 numnodes IN 总进程数, * 必须是 2^n, n = 1,2,3,...

47 IV b IN 下次对角线 a IN 对角线 c IN 上次对角线 x OUT 返回方程的解 q IN 线性方程组的右端项 */ void Thomas_P(int mynode, int numnodes, int N, double *b, double *a, double *c, double *x, double *q){ int i,j,k,i_global; int rows_local,local_offset; double S[2][2],T[2][2],s1tmp,s2tmp; double *l,*d,*y; MPI_Status status;

48 V MPI_Op op; MPI_Op_create(&oper, 0, &op); l = new double[n]; d = new double[n]; y = new double[n]; for(i=0;i<n;i++) l[i] = d[i] = y[i] = 0.0; S[0][0] = S[1][1] = 1.0; S[1][0] = S[0][1] = 0.0; // 每个处理器处理的行数 rows_local = (int) floor(n/numnodes); // 每个处理器处理的行数的第一行的位置

49 VI local_offset = mynode*rows_local; 计算每个处理器上 R k 的乘积 R k = ak b k c k // Form local products of R_k matrices if(mynode==0){ s1tmp = a[local_offset]*s[0][0]; S[1][0] = S[0][0]; S[1][1] = S[0][1]; S[0][1] = a[local_offset]*s[0][1]; S[0][0] = s1tmp; for(i=1;i<rows_local;i++){ s1tmp = a[i+local_offset]*s[0][0] -

50 VII b[i+local_offset-1]* c[i+local_offset-1]*s[1][0]; s2tmp = a[i+local_offset]*s[0][1] - b[i+local_offset-1]* c[i+local_offset-1]*s[1][1]; S[1][0] = S[0][0]; S[1][1] = S[0][1]; S[0][0] = s1tmp; S[0][1] = s2tmp; } } else{ for(i=0;i<rows_local;i++){ s1tmp = a[i+local_offset]*s[0][0] - b[i+local_offset-1]* c[i+local_offset-1]*s[1][0];

51 VIII } } s2tmp = a[i+local_offset]*s[0][1] - b[i+local_offset-1]* c[i+local_offset-1]*s[1][1]; S[1][0] = S[0][0]; S[1][1] = S[0][1]; S[0][0] = s1tmp; S[0][1] = s2tmp; 前缀归约的两种算法算法 1 (MPI_Scan ) 和算法 2

52 IX MPI_Scan(S, S, 4, MPI_DOUBLE, op, MPI_COMM_WORLD); MPI_Op_free(&op); // Full-recursive doubling algorithm // for distribution /* 这一部分被上面的 MPI_Scan 函数代替 for(i=0; i<=log2(numnodes);i++){ if(mynode+pow(2,i) < numnodes) MPI_Send(S,4,MPI_DOUBLE, int(mynode+pow(2,i)),0, MPI_COMM_WORLD); if(mynode-pow(2,i)>=0){ MPI_Recv(T,4,MPI_DOUBLE, int(mynode-pow(2,i)),0,

53 X } */ } MPI_COMM_WORLD,&status); s1tmp = S[0][0]*T[0][0] + S[0][1]*T[1][0]; S[0][1] = S[0][0]*T[0][1] + S[0][1]*T[1][1]; S[0][0] = s1tmp; s1tmp = S[1][0]*T[0][0] + S[1][1]*T[1][0]; S[1][1] = S[1][0]*T[0][1] + S[1][1]*T[1][1]; S[1][0] = s1tmp;

54 XI 计算传送 d k 把每个处理器上最后的 d k 算出来, 然后送到下一个处理器上这一段程序可以用 MPI_SendRecv 实现吗? //Calculate last d_k first so that //it can be distributed, //and then do the distribution. d[local_offset+rows_local-1] = (S[0][0] + S[0][1])/ (S[1][0] + S[1][1]); if(mynode == 0){ MPI_Send(&d[local_offset+rows_local-1],1, MPI_DOUBLE, 1,0,MPI_COMM_WORLD); }

55 XII else{ MPI_Recv(&d[local_offset-1],1,MPI_DOUBLE, mynode-1,0,mpi_comm_world,&status); if(mynode!= numnodes-1) MPI_Send(&d[local_offset+rows_local-1], 1,MPI_DOUBLE, mynode+1,0,mpi_comm_world); } 计算所有的 l k 和 d k 利用上一个处理器送来的 d k 计算本处理器上的所有 l k 和 d k

56 XIII // Compute in parallel the local values of d_k and l_k if(mynode == 0){ l[0] = 0; d[0] = a[0]; for(i=1;i<rows_local-1;i++){ l[local_offset+i] = b[local_offset+i-1]/ d[local_offset+i-1]; d[local_offset+i] = a[local_offset+i] - l[local_offset+i]*c[local_offset+i-1]; } l[local_offset+rows_local-1] = b[local_offset+rows_local-2]/ d[local_offset+rows_local-2]; } else{ for(i=0;i<rows_local-1;i++){ l[local_offset+i] = b[local_offset+i-1]/ d[local_offset+i-1]; d[local_offset+i] = a[local_offset+i] - l[local_offset+i]*c[local_offset+i-1]; } l[local_offset+rows_local-1] = b[local_offset+rows_local-2]/ d[local_offset+rows_local-2]; }

57 XIV Forward substitution and backward substition 求解下三角和上三角方程组, 这里没有使用并行, 而是把分解后的 L 和 U 传回到第一个 CPU, 然后只用第一个 CPU 进行计算下面这段程序可以并行吗? 一定需要把所有的 l k and d k 传回第一个 CPU 吗? /************************************************ if(mynode>0) d[local_offset-1] = 0; // Distribute d_k and l_k to all processes double * tmp = new double[n]; for(i=0;i<n;i++) tmp[i] = d[i]; MPI_Allreduce(tmp,d,N,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);

58 XV for(i=0;i<n;i++) tmp[i] = l[i]; MPI_Allreduce(tmp,l,N,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); delete[] tmp; if(mynode ==0){ /* Forward Substitution [L][y] = [q] */ y[0] = q[0]; for(i=1;i<n;i++) y[i] = q[i] - l[i]*y[i-1]; /* Backward Substitution [U][x] = [y] */ x[n-1] = y[n-1]/d[n-1]; for(i=n-2;i>=0;i--) x[i] = (y[i] - c[i]*x[i+1])/d[i];

59 XVI } delete[] l; delete[] y; delete[] d; return; } 调用上面程序的主程序 int main(int argc, char * argv[]){ MPI_Init(&argc, &argv); int mynode, numnodes; int N = 4; MPI_Comm_rank(MPI_COMM_WORLD, &mynode); MPI_Comm_size(MPI_COMM_WORLD, &numnodes);

60 XVII double *b = new double[n]; double *a = new double[n]; double *c = new double[n]; double *x = new double[n]; double *q = new double[n]; for(int i=0;i<n;i++){ b[i] = 1.0; a[i] = 2.0; c[i] = 0.0; q[i] = i+2; } Thomas_P(mynode, numnodes, N, b, a, c, x, q);

61 XVIII } if(mynode==0){ for(int i=0; i<n; i++){ cout << x[i] <<endl; } } MPI_Finalize(); return 1;

62 框架 1 2 Method of Undetermined Coecients BVP 3 Second-Order Initial Value Problems

63 在上面的并行算法, 我们注意到有这样一种情况 n 个进程, 进程 P i 需要向进程 P i+1 发送数据第一个进程仅仅给第二个进程发送数据, 最后一个进程仅仅接收数据, 这种情况可以用 MPI_Sendrecv 函数来实现吗? 如果能请给出一个例子研究上面给出的并行 Thomas 算法的代码, 写出简洁高效率的代码要求求解 Ly = q 和 Ux = y 的过程并行处理使用本讲义中给出的四阶隐式格式离散一维 Poisson 方程, 验证数值收敛的阶数

64 框架 1 2 Method of Undetermined Coecients BVP 3 Second-Order Initial Value Problems

65 Parallel Scientic Computing in C++ and MPI, Chapter 6.

全唐诗50

全唐诗50 ... 1... 1... 2... 2... 3... 3... 3... 4... 4... 5... 5... 6... 6... 6... 7... 7... 7... 8... 8... 8... 9 I II... 9...10...10...10...11...11...11...12...12...12...13...14...14...15...15...16...16...16...17,...17...18...18...19...19...19