算法导论第十九章 并行算法：解锁计算新维度
第十九章 并行算法：解锁计算新维度
“并行计算不是未来，而是现在。” —— David Patterson
在单核性能增长放缓的时代，并行算法成为突破计算极限的关键。本章将带你探索多核处理器、分布式系统和GPU加速的奇妙世界，揭示如何通过协同计算解决传统串行算法难以企及的大规模问题。
19.1 并行计算革命
19.1.1 从单核到并行的演变
摩尔定律的变迁：
- 1986-2004：时钟频率每年增长40%
- 2004至今：核心数量每18-24个月翻倍
19.1.2 并行计算模型对比
| 模型 | 代表技术 | 通信方式 | 适用场景 |
|---|---|---|---|
| 共享内存 | OpenMP | 内存共享 | 单机多核 |
| 消息传递 | MPI | 网络通信 | 超级计算机 |
| SIMD | GPU | 数据并行 | 图形/科学计算 |
| 数据流 | TensorFlow | 数据依赖 | 机器学习 |
19.2 并行算法设计模式
19.2.1 分治并行模式
19.2.2 流水线并行模式
// 图像处理流水线示例
#pragma omp parallel sections
{#pragma omp sectionstage1_process(image); // 阶段1：去噪#pragma omp sectionstage2_process(image); // 阶段2：边缘检测#pragma omp sectionstage3_process(image); // 阶段3：特征提取
}
19.2.3 MapReduce模式
// 词频统计MapReduce伪代码
void map(string document) {for each word in document:emit(word, 1);
}void reduce(string word, list<int> counts) {int sum = 0;for each count in counts:sum += count;emit(word, sum);
}
19.3 共享内存并行：OpenMP
19.3.1 OpenMP核心指令
#include <omp.h>
#include <stdio.h>int main() {int n = 1000000;double sum = 0.0;// 并行区域#pragma omp parallel for reduction(+:sum)for (int i = 0; i < n; i++) {double term = (i % 2 == 0) ? 1.0 : -1.0;term /= 2 * i + 1;sum += term;}double pi = 4 * sum;printf("π ≈ %.15f\n", pi);return 0;
}
19.3.2 性能优化技巧
-
循环调度策略：
schedule(static)：均匀划分schedule(dynamic)：动态分配schedule(guided)：大块优先
-
数据局部性优化：
// 缓存友好访问
#pragma omp parallel for collapse(2)
for (int i = 0; i < N; i++) {for (int j = 0; j < M; j++) {// 连续内存访问A[i][j] = B[i][j] * C[j][i]; }
}
19.4 分布式并行：MPI
19.4.1 MPI通信模式
19.4.2 矩阵乘法MPI实现
#include <mpi.h>
#include <stdio.h>#define N 1024 // 矩阵大小int main(int argc, char** argv) {int rank, size;MPI_Init(&argc, &argv);MPI_Comm_rank(MPI_COMM_WORLD, &rank);MPI_Comm_size(MPI_COMM_WORLD, &size);// 计算每个进程负责的行块int rows_per_proc = N / size;double *A_local = malloc(rows_per_proc * N * sizeof(double));double *B = malloc(N * N * sizeof(double));double *C_local = malloc(rows_per_proc * N * sizeof(double));// 主进程初始化矩阵if (rank == 0) {// 初始化A和B矩阵...}// 分发数据MPI_Scatter(A, rows_per_proc * N, MPI_DOUBLE, A_local, rows_per_proc * N, MPI_DOUBLE, 0, MPI_COMM_WORLD);MPI_Bcast(B, N * N, MPI_DOUBLE, 0, MPI_COMM_WORLD);// 本地矩阵乘法for (int i = 0; i < rows_per_proc; i++) {for (int j = 0; j < N; j++) {C_local[i * N + j] = 0;for (int k = 0; k < N; k++) {C_local[i * N + j] += A_local[i * N + k] * B[k * N + j];}}}// 收集结果MPI_Gather(C_local, rows_per_proc * N, MPI_DOUBLE, C, rows_per_proc * N, MPI_DOUBLE, 0, MPI_COMM_WORLD);// 清理资源free(A_local);free(B);free(C_local);MPI_Finalize();return 0;
}
19.4.3 MPI通信优化
| 通信模式 | 函数 | 适用场景 |
|---|---|---|
| 点对点 | MPI_Send/Recv | 不规则数据交换 |
| 集合通信 | MPI_Bcast | 广播数据 |
MPI_Reduce | 全局归约 | |
MPI_Alltoall | 全交换 | |
| 非阻塞通信 | MPI_Isend/Irecv | 重叠计算和通信 |
19.5 GPU并行计算：CUDA
19.5.1 CUDA编程模型
19.5.2 CUDA向量加法
#include <stdio.h>
#include <cuda_runtime.h>// CUDA核函数
__global__ void vectorAdd(int *a, int *b, int *c, int n) {int tid = blockIdx.x * blockDim.x + threadIdx.x;if (tid < n) {c[tid] = a[tid] + b[tid];}
}int main() {int n = 1 << 20; // 1M元素size_t size = n * sizeof(int);// 主机内存int *h_a = malloc(size);int *h_b = malloc(size);int *h_c = malloc(size);// 设备内存int *d_a, *d_b, *d_c;cudaMalloc(&d_a, size);cudaMalloc(&d_b, size);cudaMalloc(&d_c, size);// 初始化数据for (int i = 0; i < n; i++) {h_a[i] = i;h_b[i] = i * 2;}// 数据复制到设备cudaMemcpy(d_a, h_a, size, cudaMemcpyHostToDevice);cudaMemcpy(d_b, h_b, size, cudaMemcpyHostToDevice);// 启动核函数int threadsPerBlock = 256;int blocksPerGrid = (n + threadsPerBlock - 1) / threadsPerBlock;vectorAdd<<<blocksPerGrid, threadsPerBlock>>>(d_a, d_b, d_c, n);// 结果复制回主机cudaMemcpy(h_c, d_c, size, cudaMemcpyDeviceToHost);// 验证结果for (int i = 0; i < n; i++) {if (h_c[i] != h_a[i] + h_b[i]) {printf("Error at index %d\n", i);break;}}// 清理资源free(h_a); free(h_b); free(h_c);cudaFree(d_a); cudaFree(d_b); cudaFree(d_c);return 0;
}
19.5.3 CUDA性能优化
-
内存层次优化：
- 全局内存 → 共享内存 → 寄存器
- 合并内存访问
-
执行配置优化：
dim3 blocks(256, 1, 1);
dim3 grid((width+255)/256, (height+63)/64, 1);
kernel<<<grid, blocks>>>(...);
- 异步执行：
cudaMemcpyAsync(d_data, h_data, size, cudaMemcpyHostToDevice, stream);
kernel<<<grid, block, 0, stream>>>(...);
cudaMemcpyAsync(h_result, d_result, size, cudaMemcpyDeviceToHost, stream);
19.6 并行图算法
19.6.1 并行BFS实现
void parallel_bfs(Graph *g, int start) {int *visited = calloc(g->nvertices, sizeof(int));int *current_frontier = malloc(g->nvertices * sizeof(int));int *next_frontier = malloc(g->nvertices * sizeof(int));int current_size = 1;current_frontier[0] = start;visited[start] = 1;while (current_size > 0) {int next_size = 0;#pragma omp parallel for schedule(dynamic, 64)for (int i = 0; i < current_size; i++) {int node = current_frontier[i];for (int j = 0; j < g->degree[node]; j++) {int neighbor = g->edges[node][j];if (!__sync_fetch_and_or(&visited[neighbor], 1)) {int pos;#pragma omp atomic capturepos = next_size++;next_frontier[pos] = neighbor;}}}// 交换边界int *temp = current_frontier;current_frontier = next_frontier;next_frontier = temp;current_size = next_size;}free(visited);free(current_frontier);free(next_frontier);
}
19.6.2 并行PageRank算法
void parallel_pagerank(Graph *g, double damping, double epsilon) {int n = g->nvertices;double *rank = malloc(n * sizeof(double));double *new_rank = malloc(n * sizeof(double));// 初始化#pragma omp parallel forfor (int i = 0; i < n; i++) {rank[i] = 1.0 / n;}double diff;int iter = 0;do {diff = 0;#pragma omp parallel forfor (int i = 0; i < n; i++) {new_rank[i] = (1 - damping) / n;}// 传播排名#pragma omp parallel for reduction(+:diff)for (int i = 0; i < n; i++) {if (g->degree[i] > 0) {double contribution = damping * rank[i] / g->degree[i];for (int j = 0; j < g->degree[i]; j++) {int neighbor = g->edges[i][j];#pragma omp atomicnew_rank[neighbor] += contribution;}}}// 计算差异#pragma omp parallel for reduction(+:diff)for (int i = 0; i < n; i++) {diff += fabs(new_rank[i] - rank[i]);rank[i] = new_rank[i];}printf("迭代 %d, 差异: %f\n", ++iter, diff);} while (diff > epsilon);free(rank);free(new_rank);
}
19.7 并行排序算法
19.7.1 并行快速排序
void parallel_quicksort(int *arr, int left, int right) {if (right - left < 1000) {qsort(arr + left, right - left + 1, sizeof(int), compare);return;}int pivot = partition(arr, left, right);#pragma omp task shared(arr)parallel_quicksort(arr, left, pivot - 1);#pragma omp task shared(arr)parallel_quicksort(arr, pivot + 1, right);#pragma omp taskwait
}// 调用方式
#pragma omp parallel
{#pragma omp singleparallel_quicksort(arr, 0, n-1);
}
19.7.2 并行归并排序
void parallel_mergesort(int *arr, int l, int r) {if (l < r) {int m = l + (r - l) / 2;#pragma omp task if(r-l > 100000)parallel_mergesort(arr, l, m);#pragma omp task if(r-l > 100000)parallel_mergesort(arr, m+1, r);#pragma omp taskwaitmerge(arr, l, m, r);}
}
19.7.3 排序算法性能对比
| 算法 | 时间复杂度 | 并行效率 | 适用数据规模 |
|---|---|---|---|
| 快速排序 | O(n log n) | 高 | 中小规模 |
| 归并排序 | O(n log n) | 高 | 大规模 |
| 基数排序 | O(nk) | 中等 | 整数数据 |
| 样本排序 | O(n log n) | 高 | 超大规模 |
19.8 并行性能分析
19.8.1 加速比定律
-
Amdahl定律：
S = 1 ( 1 − P ) + P N S = \frac{1}{(1 - P) + \frac{P}{N}} S=(1−P)+NP1 -
Gustafson定律：
S = N + ( 1 − N ) ( 1 − P ) S = N + (1 - N)(1 - P) S=N+(1−N)(1−P)
19.8.2 性能分析工具
- 时间线分析：
$ nsys profile ./parallel_program
- 性能计数器：
#include <papi.h>int events[2] = {PAPI_TOT_CYC, PAPI_TOT_INS};
long long values[2];PAPI_start_counters(events, 2);// 并行区域...PAPI_stop_counters(values, 2);
double cpi = (double)values[0] / values[1];
19.9 实际应用案例
19.9.1 气候模拟（MPI）
void climate_simulation() {// 初始化MPIMPI_Init(...);// 划分全球网格int local_rows = global_rows / size;Grid *local_grid = create_grid(local_rows, COLS);while (time < MAX_TIME) {// 边界交换exchange_boundaries(local_grid);// 本地物理计算compute_physics(local_grid);// 全局同步MPI_Barrier(MPI_COMM_WORLD);time += TIME_STEP;}MPI_Finalize();
}
19.9.2 深度学习训练（GPU）
void train_neural_network() {// 初始化模型和数据Model *model = create_model();DataBatch *batch = next_batch();while (!converged) {// 前向传播forward_propagate<<<blocks, threads>>>(model, batch);// 损失计算compute_loss<<<blocks, threads>>>(model, batch);// 反向传播backward_propagate<<<blocks, threads>>>(model, batch);// 参数更新update_parameters<<<blocks, threads>>>(model);batch = next_batch();}
}
19.9.3 实时渲染（多核+GPU）
void render_frame() {// CPU任务：场景管理#pragma omp parallel sections{#pragma omp sectionupdate_physics();#pragma omp sectionprocess_ai();#pragma omp sectionhandle_input();}// GPU任务：渲染render_scene<<<grid, block>>>(scene);// 显示结果display_frame();
}
19.10 并行计算挑战与未来
19.10.1 当前挑战
- 负载均衡：动态任务分配
- 通信开销：减少数据传输
- 并发错误：数据竞争、死锁
- 编程复杂性：学习曲线陡峭
19.10.2 未来趋势
- 异构计算：CPU+GPU+FPGA协同
- 量子并行：量子算法加速
- 神经形态计算：模拟大脑并行处理
- 自动并行化：编译器智能优化
19.11 总结与下一章预告
并行算法将计算从线性推向多维空间：
- 模型多样：共享内存、消息传递、SIMD架构
- 框架成熟：OpenMP、MPI、CUDA各司其职
- 算法丰富：并行排序、图算法、科学计算
- 应用广泛：从科学模拟到深度学习
下一章预告：近似算法
第二十章我们将探索近似算法的智慧：
- NP难问题的实用解决方案
- 近似比与性能保证分析
- 经典问题：旅行商问题、顶点覆盖、背包问题
- 随机化近似方法
- 近似方案：PTAS与FPTAS
当精确解难以获得时，近似算法提供了在可接受时间内获得高质量解的优雅途径，平衡计算效率与结果质量。
并行计算不仅是技术演进的结果，更是解决人类面临复杂问题的关键。从天气预报到基因测序，从深度学习到宇宙模拟，并行算法正在扩展人类认知的边界，开启计算科学的新纪元。