第十九章 并行算法：解锁计算新维度
“并行计算不是未来，而是现在。” —— 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 section
stage1_process(image); // 阶段1：去噪
#pragma omp section
stage2_process(image); // 阶段2：边缘检测
#pragma omp section
stage3_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 capture
pos = 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 for
for (int i = 0; i < n; i++) {
rank[i] = 1.0 / n;
}
double diff;
int iter = 0;
do {
diff = 0;
#pragma omp parallel for
for (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 atomic
new_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 single
parallel_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 taskwait
merge(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)+PNS = \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() {
// 初始化MPI
MPI_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 section
update_physics();
#pragma omp section
process_ai();
#pragma omp section
handle_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
当精确解难以获得时，近似算法提供了在可接受时间内获得高质量解的优雅途径，平衡计算效率与结果质量。
并行计算不仅是技术演进的结果，更是解决人类面临复杂问题的关键。从天气预报到基因测序，从深度学习到宇宙模拟，并行算法正在扩展人类认知的边界，开启计算科学的新纪元。
扫一扫
6060
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
