超级会员免费看
并行计算中的多种算法与工具应用
1. MPI三维FFT示例
在高维情况下，正交方向的独立性使得并行化在概念上变得清晰。此时，不再是对块进行转置，而是移动长方体（铅笔形状）。三维变换分三个阶段计算（对应每个 $x_1$、$x_2$、$x_3$ 方向），通过对另外两个方向进行索引来实现并行化。
1.1 转置函数
以下是用于立方FFT的3 - D转置函数 Xpose 的代码：
void Xpose(float *a, int nz, int nx) {
/* 3-D transpose for cubic FFT: nx = the thickness
of the slabs, nz=dimension of the transform */
float t0,t1;
static float *buf_io;
int i, ijk, j, js, k, step, n2, nb, np, off;
static int init=-1;
int size, rk, other;
MPI_Status stat;
MPI_Comm_size(MPI_COMM_WORLD, &size);
MPI_Comm_rank(MPI_COMM_WORLD, &rk);
/* number of local planes of 3D array */
n2 = 2*nz; np = 2*nx*nx; nb = nz*n2*nx;
if(init!=nx){
if(init>0) free(buf_io);
buf_io = (float *)malloc(nb*sizeof(float));
init = nx;
}
/* local transpose of first block (in-place) */
for(j = 0; j < nx; j++){
off = j*2*nx + rk*n2;
for(k = 0; k < nz; k ++){
for(i = 0; i < k; i++) {
t0 = a[off + i*np + k*2];
t1 = a[off + i*np + k*2+1];
a[off+i*np+k*2] = a[off+k*np+2*i];
a[off+i*np+k*2+1] = a[off+k*np+2*i+1];
a[off+k*np+2*i] = t0;
a[off+k*np+2*i+1] = t1;
}
}
}
/* size-1 communication steps */
for (step = 1; step < size; step ++) {
other = rk - step;
/* fill send buffer */
ijk = 0;
for(j=0;j<nx;j++){
for(k=0;k<nz;k++){
off = j*2*nx + other*n2 + k*np;
for(i=0;i<n2;i++){
buf_io[ijk++] = a[off + i];
}
}
}
/* exchange data */
MPI_Sendrecv_replace(buf_io,n2*nz*nx,MPI_FLOAT,
other , rk , other , other , MPI_COMM_WORLD , &stat );
/* write back recv buffer in transposed order */
ijk = 0;
for(j=0;j<nx;j++){
off = j*2*nx + other*n2;
for (k=0 ; k<nz ; k++) {
for(i=0;i<nz;i++){
a[off+i*np+2*k] = buf_io[ijk];
a[off+i*np+2*k+1] = buf_io[ijk+1];
ijk += 2;
}
}
}
}
}
1.2 三维变换函数
假设数组 w 已用 $n/2$ 次单位根的幂初始化，以下是三维变换函数 FFT3D 的代码：
void FFT3D(float *a,float *w,float sign,int nz,int nx)
{
int i,j,k,off,rk;
static int nfirst=-1;
static float *pw;
float *pa;
void Xpose();
void cfft2();
MPI_Comm_rank(MPI_COMM_WORLD,&rk);
if (nfirst!=nx) {
if(nfirst>0) free(pw);
pw = (float *) malloc(2*nx*sizeof(float));
nfirst = nx;
}
/* X-direction */
for(k=0;k<nz;k++){
off = 2*k*nx*nx;
for(j=0;j<nx;j++){
pa = a + off + 2*nx*j;
cfft2(nx,pa,w,sign);
}
}
/* Y-direction */
for(k=0;k<nz;k++){
for(i=0;i<nx;i++){
off = 2*k*nx*nx+2*i;
for(j=0;j<nx;j++){
*(pw+2*j) = *(a+2*j*nx+off);
*(pw+2*j+1) = *(a+2*j*nx+1+off);
}
cfft2(nx,pw,w,sign);
for(j=0;j<nx;j++){
*(a+2*j*nx+off) = *(pw+2*j);
*(a+2*j*nx+1+off) = *(pw+2*j+1);
}
}
}
/* Z-direction */
Xpose(a,nz,nx);
for(k=0;k<nz;k++){
off = 2*k*nx*nx;
for(j=0;j<nx;j++){
pa = a + off + 2*nx*j;
cfft2(nx,pa,w,sign);
}
}
Xpose(a,nz,nx);
}
2. MPI蒙特卡罗（MC）积分示例
在MC模拟中，将工作分配到多个CPU的最简单方法是将样本分成 $p = N_{CPUs}$ 份，每个CPU计算大小为 $N/p$ 的样本。另一种方法是进行域分解，将积分域的各个部分分配给每个处理器。
2.1 代码示例
以下是对一个任意函数 $f(x, y)$ 在星形区域上进行积分的代码：
#include <stdio.h>
#include <math.h>
#include "mpi.h"
#define NP 10000
/* f(x,y) is any "reasonable" function */
#define f(x,y) exp(-x*x-0.5*y*y)*(10.0*x*x-x)*(y*y-y)
main(int argc, char **argv)
{
/* Computes integral of f() over star-shaped
region: central square is 1 x 1, each of
N,S,E,W points are 2 x 1 triangles of same
area as center square */
int i,ierr,j,size,ip,master,rank;
/* seeds for each CPU: */
float seeds[]={331.0,557.0,907.0,1103.0,1303.0};
float x,y,t,tot;
static float seed;
float buff[1];
/* buff for partial sums */
float sum[4];
/* part sums from "points" */
float ggl();
/* random number generator */
MPI_Status stat;
MPI_Init(&argc,&argv);
/* initialize MPI */
MPI_Comm_size(MPI_COMM_WORLD,&size); /* 5 cpus */
MPI_Comm_rank(MPI_COMM_WORLD,&rank); /* rank */
master =0;
/* master */
ip = rank;
if (ip==master){
/* master computes integral over center square */
tot = 0.0;
seed = seeds[ip];
for(i=0;i<NP;i++){
x = ggl(&seed)-0.5;
y = ggl(&seed)-0.5;
tot += f(x,y);
}
tot *= 1.0/((float) NP); /* center part */
} else {
/* rank != 0 computes E,N,W,S points of star */
seed = seeds [ip];
sum[ip-1] = 0.0;
for(i=0;i<NP;i++){
x = ggl(&seed);
y = ggl(&seed) - 0.5;
if(y > (0.5-0.25*x)){
x = 2.0 - x; y = -(y-0.5);
}
if(y < (-0.5+0.25*x)){
x = 2.0 - x; y = -(y+0.5);
}
x += 0.5;
if(ip==2){
t = x; x = y; y = -t;
} else if(ip==3){
x = -x; y = -y;
} else if(ip==4){
t = x;
x = -y; y = t;
}
sum[ip-1] += f(x,y);
}
sum[ip-1] *= 1.0/((float) NP);
buff[0] = sum[ip-1];
MPI_Send(buff,1,MPI_FLOAT,0,0,MPI_COMM_WORLD);
}
if(ip==master){ /* get part sums of other cpus */
for(i=1;i<5;i++){
MPI_Recv(buff,1,MPI_FLOAT,MPI_ANY_SOURCE,
MPI_ANY_TAG,MPI_COMM_WORLD,
&stat);
tot += buff[0];
}
printf(" Integral = %e\n",tot);
}
MPI_Finalize();
}
3. PETSc简介
PETSc是用于科学计算的可移植扩展工具包，是一组用于在并行（和串行）计算机上求解大型稀疏线性和非线性方程组的C例程。它使用MPI标准进行所有消息传递通信。
3.1 矩阵和向量
3.1.1 矩阵创建
矩阵通过 MatCreate 命令创建：
Mat A;
/* linear system matrix */
ierr = MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,
PETSC_DECIDE,n*n,n*n,&A);
ierr = MatSetFromOptions(A);
/*
Set up the system matrix */
ierr = MatGetOwnershipRange(A,&Istart,&Iend);
for (I=Istart; I<Iend; I++) {
v = -1.0; i = I/n; j = I - i*n;
if (i>0) {J = I - n;
ierr = MatSetValues(A,1,&I,1,&J,&v,
INSERT_VALUES);
}
v = 4.0;
ierr = MatSetValues(A, 1 ,&I, 1,&I,&v,
INSERT_VALUES);
}
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
3.1.2 向量创建
向量的创建方式与矩阵类似：
Vec u;
ierr = VecCreate(PETSC_COMM_WORLD,&u);
ierr = VecSetSizes(u,PETSC_DECIDE,n*n);
ierr = VecSetFromOptions(u);
3.2 Krylov子空间方法和预条件器
3.2.1 线性求解器上下文创建
/* Create linear solver context */
ierr = SLESCreate(PETSC_COMM_WORLD,&sles);
/*
Set operator */
ierr = SLESSetOperators(sles,A,A,
DIFFERENT_NONZERO_PATTERN);
/*
Set Krylov subspace method
*/
ierr = SLESGetKSP(sles,&ksp);
ierr = KSPSetType(ksp,KSPCG);
ierr = KSPSetTolerances(ksp,1.e-6,PETSC_DEFAULT,
PETSC_DEFAULT,PETSC_DEFAULT);
3.2.2 预条件器定义
ierr = SLESGetPC(sles,&pc);
ierr = PCSetType(pc,PCJACOBI);
3.2.3 调用PETSc求解器
ierr = SLESSolve(sles,b,x,&its);
4. 数值实验
4.1 实验设置
在Beowulf PC集群上使用5 - 点有限差分方案求解泊松方程 $-\Delta u = f$。实验针对三种问题规模进行，分别对应 $n = 127$、$255$ 和 $511$。
4.2 实验结果
| n | P | Jacobi | Block Jacobi (1) | Block Jacobi (2) | IC(0) |
|---|---|---|---|---|---|
| 127 | 1 | 2.95 (217) | 4.25 (168) | 0.11 (1) | 5.29 (101) |
| 127 | 2 | 2.07 (217) | 2.76 (168) | 0.69 (22) | 4.68 (147) |
| 127 | 4 | 1.38 (217) | 1.66 (168) | 0.61 (38) | 3.32 (139) |
| 127 | 8 | 1.32 (217) | 1.53 (168) | 0.49 (30) | 3.12 (137) |
| 127 | 16 | 1.51 (217) | 1.09 (168) | 0.50 (65) | 2.02 (128) |
| 255 | 1 | 29.0 (426) | 34.6 (284) | 0.50 (1) | 42.0 (197) |
| 255 | 2 | 17.0 (426) | 22.3 (284) | 4.01 (29) | 31.6 (263) |
| 255 | 4 | 10.8 (426) | 12.5 (284) | 3.05 (46) | 20.5 (258) |
| 255 | 8 | 6.55 (426) | 6.91(284) | 2.34 (66) | 14.6 (243) |
| 255 | 16 | 4.81 (426) | 4.12(284) | 16.0 (82) | 10.9 (245) |
| 255 | 32 | 4.23 (426) | 80.9 (284) | 2.09 (113) | 21.7 (241) |
| 511 | 1 | 230.0 (836) | 244.9 (547) | 4.1 (1) | 320.6 (384) |
| 511 | 2 | 152.2 (836) | 157.3 (547) | 36.2 (43) | 253.1 (517) |
| 511 | 4 | 87.0 (836) | 86.2 (547) | 25.7 (64) | 127.4 (480) |
| 511 | 8 | 54.1 (836) | 53.3 (547) | 15.1 (85) | 66.5 (436) |
| 511 | 16 | 24.5 (836) | 24.1 (547) | 21.9 (110) | 36.6 (422) |
| 511 | 32 | 256.8 (836) | 17.7 (547) | 34.5 (135) | 107.9 (427) |
4.3 结果分析
- 前两列的预条件器不依赖于处理器数量，因此对于给定问题规模，迭代次数是恒定的。
- 在
Block Jacobi (2)和IC(0)情况下，预条件器依赖于处理器数量，迭代次数也会随 $p$ 变化。 - 对于小问题规模，加速效果不明显。对于较大问题，最多使用16个处理器仍能获得不断增加的加速比。
- 在Beowulf集群上，
Block Jacobi (2)预条件器效果最佳，它考虑了所有局部信息并进行直接求解，迭代次数最少。IC(0)虽然所需内存较少，但在减少迭代次数方面不如Block Jacobi (2)。
5. 练习题
5.1 有效带宽测试
从 Pallas 下载 EFF_BW 基准测试。具体步骤如下：
1. 解压Pallas EFF_BW 测试。
2. 编辑机器的 makefile 并构建基准测试。
3. 针对不同的消息长度和CPU数量运行测试，确定：
- 各种消息大小的带宽。
- 仅进行握手的延迟，即外推到零消息长度以获得延迟。
5.2 并行MC积分
修改并行MC数值积分代码，具体操作如下：
1. 测试独立随机数流的重要性：
- 更改代码中每个CPU的随机数种子，观察对积分结果的影响。
- 在每个处理器上使用相同的种子，观察积分结果是否有明显变化。
- 总结独立流对数值积分的重要性。
2. 进一步细分每个独立区域，修改代码以在20个处理器上运行。
3. 在20个CPU上运行积分并与5 - CPU版本的结果进行比较，再次测试对随机数流的依赖性。
5.3 用MC求解偏微分方程：第二部分
使用MPI并行化之前求解椭圆偏微分方程的MC方法，具体步骤如下：
1. 从之前的解决方案开始，并行计算多个初始 $x$ 值。
2. 可以通过编写PBS脚本或进行思想实验，想象多个批处理作业并行计算不同的 $w(x)$ 值。
5.4 用PBLAS进行矩阵 - 向量乘法
初始化分布式数据（矩阵和两个向量）并使用并行BLAS的一个例程，具体步骤如下：
1. 初始化一个尽可能方形的 $p_r \times p_c$ 进程网格，使得 $p_r \cdot p_c = p$（所需处理器数量）。
2. 将一个 $50 \times 100$ 的矩阵 $A$ 以 $5 \times 5$ 块的循环块方式分布在进程网格上，并初始化矩阵。
3. 将向量 $x$ 分布在进程网格的第一行并初始化，使得 $x_i = (-1)^i$。
4. 将结果存储在分布在进程网格第一列的向量 $y$ 中。
5. 调用PBLAS矩阵 - 向量乘法 pdgemv 计算 $y = Ax$。
5.5 使用ScaLAPACK求解 $Ax = b$
继续上一个练习，求解方程组 $Ax = b$，具体步骤如下：
1. 确保矩阵 $A$ 是方阵且非奇异，向量 $x$ 和 $y$ 长度相同。
2. 随机初始化矩阵 $A$ 的元素，将对角线元素设置为 $n$，将向量 $x$ 的所有元素设置为1。
3. 计算 $y = Ax$。
4. 调用ScaLAPACK子例程 pdgesv 求解 $Ax = y$。
5. 使用适当的函数检查结果是否与初始向量 $x$ 一致。
5.6 通过主 - 从过程分布矩阵
编写一个使用MPI原语的函数，将数据从进程0（主进程）分布到所有参与进程，具体步骤如下：
1. 使用 Chapter5/uebung6.tar 作为工作的起点，其中包含骨架和 makefile 。
2. 检查 qsub_mpich 的参数列表并在C代码中使用这些参数。
3. 尽量不浪费各个进程的本地存储空间，确保ScaLAPACK所需的所有本地内存块大小相同。
5.7 更多ScaLAPACK内容
测量 pdgetrf 和 pdgetrs 的执行时间和加速比，并分析它们不同的原因。具体步骤如下：
1. 运行求解方程组的代码，记录 pdgetrf 和 pdgetrs 的执行时间。
2. 计算不同处理器数量下的加速比。
3. 分析两者执行时间和加速比不同的原因。
6. 练习题详细分析与解答
6.1 有效带宽测试
操作步骤详解
- 解压测试文件 ：从指定的Pallas链接下载
EFF_BW基准测试文件后，使用相应的解压命令（如tar -xvf对于.tar文件）将其解压到指定目录。 - 编辑
makefile：打开解压后的makefile文件，根据自己机器的环境进行必要的修改，例如编译器路径、库文件路径等。修改完成后，在终端中执行make命令来构建基准测试程序。 - 运行测试 ：编写脚本或使用命令行工具，针对不同的消息长度和CPU数量运行测试程序。在每次运行时，记录下消息大小、带宽和延迟等相关数据。可以使用循环结构来自动化这个过程，例如：
for message_size in 1 2 4 8 16 32 64 128; do
for num_cpus in 1 2 4 8 16; do
./eff_bw -s $message_size -n $num_cpus
done
done
- 数据分析 ：对记录的数据进行处理，使用绘图工具（如Python的
matplotlib库）绘制带宽随消息大小变化的曲线，以及延迟随消息大小变化的曲线。通过线性外推的方法，将延迟曲线外推到零消息长度，得到仅进行握手的延迟。
6.2 并行MC积分
测试独立随机数流的重要性
- 更改随机数种子 ：打开代码文件，找到定义随机数种子的部分，将种子值修改为除零以外的其他值。例如：
float seeds[]={123.0, 456.0, 789.0, 1011.0, 1213.0};
编译并运行代码，记录积分结果。多次更改种子值并重复上述过程，观察积分结果的变化。
- 使用相同种子 ：将所有处理器的种子值设置为相同的值，例如：
float seeds[]={1.0, 1.0, 1.0, 1.0, 1.0};
编译并运行代码，与之前使用不同种子的结果进行比较，观察是否有明显差异。
- 总结重要性 ：根据实验结果，总结独立随机数流对数值积分的重要性。如果使用相同种子时积分结果与使用不同种子时有明显差异，说明独立随机数流对于准确的数值积分非常重要。
细分独立区域并在20个处理器上运行
- 分析区域细分方法 ：参考图5.24，对星形区域的每个部分进行进一步细分。例如，将中心正方形细分为四个相等的小正方形，将每个三角形细分为四个相似的小三角形。
- 修改代码 ：根据细分方法，修改代码中的采样逻辑和处理器分配逻辑。确保每个处理器负责一个细分区域的积分计算。例如，对于中心正方形的细分，可以在代码中添加额外的条件判断来分配不同的采样范围。
if (rank < 4) {
// 处理中心正方形细分区域
float x_offset = (rank % 2) * 0.5;
float y_offset = (rank / 2) * 0.5;
for (i = 0; i < NP; i++) {
x = ggl(&seed) * 0.5 + x_offset - 0.25;
y = ggl(&seed) * 0.5 + y_offset - 0.25;
tot += f(x, y);
}
} else {
// 处理星形其他部分
...
}
- 编译和运行 ：编译修改后的代码，并在20个处理器上运行。记录积分结果。
结果比较与再次测试
- 结果比较 ：将20个CPU运行的积分结果与5 - CPU版本的结果进行比较，计算两者的相对误差。如果相对误差在可接受的范围内，说明细分区域和并行计算的方法是有效的。
- 再次测试随机数流依赖性 ：重复上述更改随机数种子和使用相同种子的实验，观察20个CPU运行时积分结果对随机数流的依赖性。
6.3 用MC求解偏微分方程：第二部分
并行计算多个初始 $x$ 值
- 修改代码 ：在之前求解椭圆偏微分方程的MC方法代码基础上，使用MPI的并行机制来计算多个初始 $x$ 值。可以使用
MPI_Scatter函数将初始 $x$ 值分配给不同的处理器，每个处理器独立计算对应的 $w(x)$ 值。
#include <mpi.h>
// 假设初始x值存储在x_values数组中
float x_values[NUM_X_VALUES];
float local_x;
MPI_Scatter(x_values, 1, MPI_FLOAT, &local_x, 1, MPI_FLOAT, 0, MPI_COMM_WORLD);
// 计算w(local_x)
float w_local = compute_w(local_x);
- 收集结果 ：使用
MPI_Gather函数将各个处理器计算得到的 $w(x)$ 值收集到主处理器进行汇总。
float w_results[NUM_X_VALUES];
MPI_Gather(&w_local, 1, MPI_FLOAT, w_results, 1, MPI_FLOAT, 0, MPI_COMM_WORLD);
批处理作业并行计算
- 编写PBS脚本 ：PBS（Portable Batch System）是一种常用的作业调度系统。编写一个PBS脚本，在脚本中指定多个批处理作业，每个作业计算一个不同的 $w(x)$ 值。脚本示例如下：
#!/bin/bash
#PBS -N MC_PDE
#PBS -l nodes=1:ppn=1
#PBS -l walltime=01:00:00
cd $PBS_O_WORKDIR
./mc_pde -x 1.0
./mc_pde -x 2.0
./mc_pde -x 3.0
...
- 提交脚本 ：将编写好的PBS脚本提交到作业调度系统中，系统会自动并行执行这些批处理作业。
6.4 用PBLAS进行矩阵 - 向量乘法
初始化进程网格
- 确定进程网格大小 ：根据所需的处理器数量 $p$，找到尽可能方形的 $p_r$ 和 $p_c$，使得 $p_r \cdot p_c = p$。可以使用循环结构来实现：
int p;
MPI_Comm_size(MPI_COMM_WORLD, &p);
int pr, pc;
for (pr = 1; pr <= p; pr++) {
if (p % pr == 0) {
pc = p / pr;
if (abs(pr - pc) < abs(pr_best - pc_best)) {
pr_best = pr;
pc_best = pc;
}
}
}
- 初始化进程网格 ：使用PBLAS提供的函数初始化进程网格。
int ictxt;
blacs_get(-1, 0, &ictxt);
blacs_gridinit(&ictxt, 'R', pr_best, pc_best);
分布矩阵和向量
- 分布矩阵 ：将 $50 \times 100$ 的矩阵 $A$ 以 $5 \times 5$ 块的循环块方式分布在进程网格上。使用PBLAS的函数进行矩阵的初始化和分布。
int descA[9];
pdgeinit(descA, ictxt, 50, 100, 5, 5, 0, 0, 50, 50);
for (int i = 0; i < 50; i++) {
for (int j = 0; j < 100; j++) {
float value = (float)(i * 100 + j);
pdgeset(descA, i, j, value);
}
}
- 分布向量 $x$ ：将向量 $x$ 分布在进程网格的第一行并初始化，使得 $x_i = (-1)^i$。
int descX[9];
pdgeinit(descX, ictxt, 100, 1, 5, 1, 0, 0, 100, 100);
for (int i = 0; i < 100; i++) {
float value = (i % 2 == 0) ? 1.0 : -1.0;
pdgeset(descX, i, 0, value);
}
- 分布向量 $y$ ：将结果向量 $y$ 分布在进程网格的第一列。
int descY[9];
pdgeinit(descY, ictxt, 50, 1, 5, 1, 0, 0, 50, 50);
调用矩阵 - 向量乘法函数
使用PBLAS的 pdgemv 函数计算 $y = Ax$。
pdgemv('N', 50, 100, 1.0, descA, 0, 0, descX, 0, 0, 0.0, descY, 0, 0);
6.5 使用ScaLAPACK求解 $Ax = b$
矩阵和向量初始化
- 初始化矩阵 $A$ ：随机初始化矩阵 $A$ 的元素，将对角线元素设置为 $n$。
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
if (i == j) {
pdgeset(descA, i, j, (float)n);
} else {
pdgeset(descA, i, j, (float)rand() / RAND_MAX);
}
}
}
- 初始化向量 $x$ ：将向量 $x$ 的所有元素设置为1。
for (int i = 0; i < n; i++) {
pdgeset(descX, i, 0, 1.0);
}
计算 $y = Ax$
使用PBLAS的 pdgemv 函数计算 $y = Ax$。
pdgemv('N', n, n, 1.0, descA, 0, 0, descX, 0, 0, 0.0, descY, 0, 0);
求解 $Ax = y$
调用ScaLAPACK的 pdgesv 函数求解 $Ax = y$。
int ipiv[n];
pdgesv(n, 1, descA, 0, 0, ipiv, descY, 0, 0, &info);
检查结果
使用适当的函数（如计算向量的误差范数）检查求解结果是否与初始向量 $x$ 一致。
float error = 0.0;
for (int i = 0; i < n; i++) {
float diff = pdgeget(descY, i, 0) - pdgeget(descX, i, 0);
error += diff * diff;
}
error = sqrt(error);
if (error < 1e-6) {
printf("Solution is correct.\n");
} else {
printf("Solution has an error: %f\n", error);
}
6.6 通过主 - 从过程分布矩阵
读取数据
主进程（进程0）从文件中读取矩阵数据到二维数组 Aglob 和向量 bglob 中。
if (rank == 0) {
FILE *file = fopen("matrix_data.txt", "r");
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
fscanf(file, "%f", &Aglob[i][j]);
}
fscanf(file, "%f", &bglob[i]);
}
fclose(file);
}
分布数据
使用MPI的通信函数将数据从主进程分布到其他进程。可以使用 MPI_Scatterv 函数来实现。
int sendcounts[p];
int displs[p];
// 计算每个进程的发送数量和位移
for (int i = 0; i < p; i++) {
sendcounts[i] = local_size;
displs[i] = i * local_size;
}
// 分布矩阵数据
MPI_Scatterv(&Aglob[0][0], sendcounts, displs, MPI_FLOAT, &Alocal[0][0], local_size, MPI_FLOAT, 0, MPI_COMM_WORLD);
// 分布向量数据
MPI_Scatterv(bglob, sendcounts, displs, MPI_FLOAT, blocal, local_size, MPI_FLOAT, 0, MPI_COMM_WORLD);
6.7 更多ScaLAPACK内容
测量执行时间
在求解方程组的代码中，使用 MPI_Wtime 函数记录 pdgetrf 和 pdgetrs 的执行时间。
double start_time, end_time;
// 测量pdgetrf的执行时间
start_time = MPI_Wtime();
pdgetrf(n, n, descA, 0, 0, ipiv, &info);
end_time = MPI_Wtime();
double time_getrf = end_time - start_time;
// 测量pdgetrs的执行时间
start_time = MPI_Wtime();
pdgetrs('N', n, 1, descA, 0, 0, ipiv, descY, 0, 0, &info);
end_time = MPI_Wtime();
double time_getrs = end_time - start_time;
计算加速比
在不同的处理器数量下运行代码，记录每次运行的执行时间，计算加速比。加速比的计算公式为：$S = \frac{T_1}{T_p}$，其中 $T_1$ 是单处理器的执行时间，$T_p$ 是 $p$ 个处理器的执行时间。
分析差异原因
- 算法复杂度 ：
pdgetrf是矩阵的LU分解过程，其算法复杂度较高，通常为 $O(n^3)$。而pdgetrs是前向和后向替换过程，算法复杂度为 $O(n^2)$。因此，pdgetrf的执行时间通常会比pdgetrs长。 - 通信开销 ：在并行计算中，
pdgetrf可能需要更多的处理器间通信来完成矩阵的分解，而pdgetrs的通信开销相对较小。这也会导致两者的执行时间和加速比不同。
7. 总结
本文介绍了并行计算中的多种算法和工具，包括MPI三维FFT、MPI蒙特卡罗积分、PETSc工具包等。通过具体的代码示例和数值实验，展示了这些算法和工具的使用方法和性能特点。同时，详细分析了练习题的解答过程，帮助读者更好地理解和掌握并行计算的相关知识。在实际应用中，需要根据具体问题的特点选择合适的算法和工具，并进行合理的并行化设计，以提高计算效率和性能。
7.1 算法和工具的优势
- MPI三维FFT ：利用正交方向的独立性，通过转置和并行计算实现了高效的三维快速傅里叶变换。
- MPI蒙特卡罗积分 ：通过样本分割和域分解的方法，将积分计算任务分配到多个处理器上，减少了通信开销，提高了计算效率。
- PETSc ：提供了丰富的矩阵和向量操作函数，以及Krylov子空间方法和预条件器，方便用户求解大型稀疏线性和非线性方程组。
7.2 注意事项
- 随机数流 ：在并行蒙特卡罗积分中，独立的随机数流对于准确的数值积分非常重要。
- 处理器数量 ：对于不同的问题规模，需要选择合适的处理器数量以获得最佳的加速比。在小问题规模下，过多的处理器可能会导致通信开销增加，反而降低性能。
- 预条件器选择 ：不同的预条件器在不同的问题和处理器数量下表现不同，需要根据具体情况进行选择。
通过对这些算法和工具的学习和实践，读者可以更好地应对并行计算中的各种挑战，提高自己的编程能力和解决实际问题的能力。
8. 流程图示例
8.1 MPI三维FFT计算流程
graph TD;
A[开始] --> B[初始化MPI];
B --> C[进行局部转置];
C --> D[进行通信步骤];
D --> E[X方向FFT];
E --> F[Y方向FFT];
F --> G[Z方向FFT];
G --> H[结束];
8.2 并行MC积分流程
graph TD;
A[开始] --> B[初始化MPI];
B --> C{是否为主进程};
C -- 是 --> D[计算中心正方形积分];
C -- 否 --> E[计算星形其他部分积分];
D --> F[接收其他进程结果];
E --> G[发送部分积分结果];
F --> H[汇总结果并输出];
G --> F;
H --> I[结束];
8.3 PETSc求解线性方程组流程
graph TD;
A[开始] --> B[初始化PETSc];
B --> C[创建矩阵和向量];
C --> D[设置矩阵元素];
D --> E[创建线性求解器上下文];
E --> F[设置算子和Krylov子空间方法];
F --> G[设置预条件器];
G --> H[调用求解器];
H --> I[输出结果];
I --> J[结束];
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
