基于 Pthread 的 CSR 稀疏矩阵 CG 并行求解

这篇文章由《并行计算》课程实验报告整理而来，保留实验思路、实现方法、核心代码、性能数据和图表，便于后续复盘。

一、实验名称与内容

1.1 实验名称

基于 Pthread 的 CSR 稀疏矩阵共轭梯度法并行实现与性能比较

1.2 实验内容

本实验以稀疏线性方程组求解为对象，采用 CSR 稀疏矩阵存储格式，实现基于共轭梯度法（Conjugate Gradient, CG）的串行程序与多线程并行程序。实验输入包括稀疏矩阵文件 matrix.txt、右端向量文件 vector.txt 以及线程数 T，程序目标为求解线性方程组：

Ax = b

其中矩阵 A 采用 CSR 格式存储，核心计算包括：

• 稀疏矩阵向量乘 SpMV
• 向量更新 x, r, p
• 点积计算 p·Ap 和 r·r

在此基础上，实验实现了 1 个串行基线程序和 3 个并行版本：

• sparse：串行版 CG 求解程序，用于正确性验证和性能基线对照
• sparse_parallel_block_reduce：按行连续分块，线程先计算局部部分和，再统一规约
• sparse_parallel_block_lock：按行连续分块，线程通过互斥锁直接累加全局点积结果
• sparse_parallel_cyclic_reduce：按循环方式划分行，线程先计算局部部分和，再统一规约

实验要求在保证结果正确的前提下，对不同并行方案在不同线程数下的运行时间进行测试，比较其加速效果与并行效率，并分析任务划分方式和规约方式对程序性能的影响。

二、实验环境的配置参数

2.1 硬件环境

根据实验平台已有说明材料，实验运行在学院提供的高性能计算集群上，主要配置如下：

• 集群结构：1 个管理节点 + 7 个计算节点
• CPU：2 x Intel(R) Xeon(R) Gold 5218
• 内存：256 GB
• 硬盘：3 x 600 GB
• 存储方式：RAID5
• 可用存储容量：约 1.2 TB

2.2 软件与运行环境

• 编程语言：C++
• 并行库：POSIX Threads（pthread）
• 作业调度环境：PBS
• 作业提交命令：qsub
• 作业查询命令：qstat

从当前实验所用 PBS 脚本可以看出，程序均在单节点共享内存环境下运行。例如：

#PBS -l nodes=1:ppn=8

因此，本实验属于单节点共享内存并行实验，主要考察线程级并行性能，而不是多节点分布式并行性能。

2.3 实验文件与输入规模

当前实验目录中的核心文件包括：

• sparse.cpp
• sparse_parallel_block_reduce.cpp
• sparse_parallel_block_lock.cpp
• sparse_parallel_cyclic_reduce.cpp
• generate_input.cpp
• matrix.txt
• vector.txt

当前测试数据规模设置为：

• 矩阵维数：n = 5000000
• 非零元个数：nnz = 34999988
• 矩阵类型：带状稀疏对称正定矩阵

2.4 说明

现有材料能够确认实验所依赖的并行编程环境、作业调度方式以及节点级硬件配置，但未给出更完整的处理器缓存参数、内存带宽实测值等细节。因此在正式报告中，应以已有材料和实际代码、脚本为依据，不对未明确给出的平台参数作额外假设。

三、实验题目问题分析

3.1 问题性质分析

本实验研究的问题是稀疏线性方程组的数值求解。给定稀疏矩阵 A 和向量 b，程序通过共轭梯度法迭代求解向量 x。与稠密矩阵问题相比，稀疏矩阵的主要特点在于：

• 非零元数量远少于 n^2
• 适合采用压缩存储格式降低空间开销
• 计算开销主要集中在 SpMV 和向量级操作

共轭梯度法适用于对称正定矩阵，迭代过程中每一步都包含矩阵向量乘、若干向量更新和点积规约。这些操作结构清晰、重复性强，适合作为共享内存并行实验对象。

3.2 可并行性分析

从程序结构看，CG 迭代中的大部分计算都可以按向量下标或矩阵行进行划分。具体来说：

• SpMV 中每一行的结果 Ap[i] 可独立计算
• 向量更新 x[i]、r[i]、p[i] 可按下标并行处理
• 点积 p·Ap 和 r·r 可先分线程求局部和，再进行全局规约

因此，本题具有较好的数据并行特征。其并行实现的关键不在于“能否并行”，而在于：

1. 如何划分矩阵行，使各线程工作量尽量均衡
2. 如何完成点积结果的全局合并，并尽量降低同步开销

3.3 可选并行方案分析

围绕上述两个关键问题，本实验主要比较了两类设计维度。

第一类是任务划分方式：

• 连续分块划分：每个线程负责一段连续矩阵行，优点是实现简单、局部性较好
• 循环划分：线程按 tid, tid+T, tid+2T, ... 的方式处理行，优点是更容易平衡负载

第二类是规约方式：

• 局部规约后统一合并：每个线程先计算局部点积和，再由主线程统一求和
• 加锁直接累加：每个线程算出局部值后，通过互斥锁直接更新全局变量

由此形成了本实验的三个并行实现：

• block + reduce
• block + lock
• cyclic + reduce

这些方案的对比重点分别是：

• 连续分块与循环划分对缓存局部性和负载均衡的影响
• 局部规约与加锁规约对同步成本的影响

3.4 本实验关注的核心问题

本实验真正关注的并不是“并行版是否能运行”，而是以下几个问题：

• 稀疏矩阵 CG 求解在共享内存环境下能获得多大加速
• 不同任务划分方式对程序性能有何影响
• 不同规约方式的同步开销差异是否明显
• 在 1 / 2 / 4 / 8 / 16 线程下，并行效率如何变化

因此，本实验既是一个数值计算并行化实验，也是一个共享内存线程划分与规约策略的对比实验。

四、方案设计

4.1 串行算法设计

本实验的并行程序是在串行 CG 求解程序基础上扩展得到的，因此首先需要明确串行算法流程。串行版本对应 sparse.cpp，其核心由 spmv_csr 和 cg_serial 两部分组成。

其中，spmv_csr 负责按 CSR 存储格式完成稀疏矩阵向量乘；cg_serial 负责完成整个共轭梯度迭代。串行程序的基本步骤为：

1. 读入 matrix.txt 中的 n、nnz、values、row_ptr 和 col_idx
2. 读入 vector.txt 中的右端向量 b
3. 初始化 x = 0，r = b，p = r
4. 计算初始残差平方和 rho0 = r·r
5. 进入迭代，计算 Ap = A * p
6. 计算点积 p·Ap，进而求出步长 alpha
7. 更新 x 与 r
8. 计算新的残差平方和 rho_new = r·r
9. 判断是否满足收敛条件
10. 由 beta = rho_new / rho0 更新搜索方向 p
11. 重复上述步骤直到收敛或达到最大迭代次数

因此，串行版已经给出了完整的算法基线。后续并行化工作并不是改变算法公式，而是在保持串行计算逻辑不变的前提下，对其中耗时较大的阶段进行线程并行。

4.2 基于串行算法的并行化思路

对串行代码进行分析可以发现，CG 迭代中的大部分工作都具有明显的数据并行特征，主要包括：

• 稀疏矩阵向量乘 Ap = A * p
• 向量更新 x[i]、r[i]、p[i]
• 点积 p·Ap 与 r·r

其中，SpMV 和向量更新可以直接按矩阵行或向量下标划分给不同线程；而点积虽然最终需要一个全局结果，但仍可以先由各线程计算局部部分和，再进行规约。

因此，当前并行代码都采用了相同的总体思路：

1. 保留串行 CG 的整体迭代框架
2. 将每轮迭代中的行级和向量级操作分配给多个线程
3. 在线程之间共享矩阵、向量和若干全局标量
4. 在阶段之间通过 pthread_barrier_t 同步
5. 对 p·Ap 和 r·r 采用不同规约策略进行比较

从设计角度看，并行版本与串行版本在数学意义上求解的是同一个问题，只是在实现层面增加了线程划分、阶段同步和结果规约机制。

4.3 三种并行方案设计

在串行算法基础上，本实验进一步设计了三种并行实现，分别从“任务划分方式”和“规约方式”两个角度展开对比。

1. 连续分块 + 局部规约

sparse_parallel_block_reduce 采用连续分块划分任务。矩阵的 n 行被尽可能平均地分给各线程，每个线程处理一个连续区间 [start, end)。这种设计的特点是：

• 每个线程处理连续行，访存局部性较好
• Ap、x、r、p 的更新都限定在各自区间内，逻辑清晰
• 线程之间没有直接写冲突

对于全局点积量 p·Ap 和 r·r，该方案采用局部规约方式。每个线程先计算自己的局部部分和，再写入线程私有位置，最后由 0 号线程统一求和并计算 alpha、beta。这种写法避免了频繁加锁，因此通常也是性能较好的实现方式。

2. 连续分块 + 加锁规约

sparse_parallel_block_lock 在任务划分上与前一方案完全一致，同样采用连续分块。它的主要区别在于点积结果的合并方式：

• 每个线程先计算局部 local_pAp
• 每个线程先计算局部 local_rho_new
• 随后通过互斥锁对全局变量进行加锁累加

与局部规约相比，这种方法实现直接，但会引入互斥锁同步开销。因此，该方案的主要作用是作为“加锁规约”的对照组，用于分析锁竞争对程序性能的影响。

3. 循环划分 + 局部规约

sparse_parallel_cyclic_reduce 保留了局部规约策略，但改变了任务划分方式。该方案中，线程从自身编号开始，按步长 num_threads 循环访问数据，即：

• 线程 0 处理 0, T, 2T, ...
• 线程 1 处理 1, T+1, 2T+1, ...
• 其余线程以此类推

因此它与 block_reduce 的区别不在于规约，而在于：

• block_reduce 采用连续分块
• cyclic_reduce 采用循环划分

循环划分的优点是更容易把工作量均匀分散到各线程；缺点是访问模式不连续，缓存局部性通常不如连续分块。该方案的实验意义主要在于比较两种划分方式对性能的影响。

4. 三种并行方案与串行基线的关系

综合来看，本实验形成了如下对照关系：

• sparse：串行基线
• sparse_parallel_block_reduce：连续分块 + 局部规约
• sparse_parallel_block_lock：连续分块 + 加锁规约
• sparse_parallel_cyclic_reduce：循环划分 + 局部规约

因此，实验比较的逻辑为：

• 先用 sparse 提供串行时间基线
• 再比较并行程序在不同线程数下的总体加速效果
• 最后在并行程序内部比较“连续分块与循环划分”“局部规约与加锁规约”的差异

4.4 程序流程图

从代码结构看，三个并行版本都不是简单地把整个 CG 交给线程独立执行，而是把一次迭代拆成多个并行阶段，并在阶段之间通过 pthread_barrier_wait 同步。根据当前实现绘制得到的程序流程图如图 4-1 所示：

图中给出了当前并行程序的总体结构。其中：

• block_reduce 和 cyclic_reduce 采用局部数组规约
• block_lock 采用互斥锁直接累加全局量
• 三个版本都共享同一套 CG 迭代框架

4.5 伪代码设计

为便于说明程序设计思路，下面给出与当前代码对应的串行与并行伪代码。

1. 串行算法伪代码

读入 CSR 矩阵 A 和向量 b
初始化 x = 0
初始化 r = b
初始化 p = r

rho0 = dot(r, r)
如果 sqrt(rho0) < tol:
    结束

for iter = 0 .. max_iter-1:
    Ap = SpMV(A, p)

    pAp = dot(p, Ap)
    如果 |pAp| 很小:
        跳出循环

    alpha = rho0 / pAp

    for i = 0 .. n-1:
        x[i] = x[i] + alpha * p[i]
        r[i] = r[i] - alpha * Ap[i]

    rho_new = dot(r, r)
    res_norm = sqrt(rho_new)
    如果 res_norm < tol:
        跳出循环

    beta = rho_new / rho0

    for i = 0 .. n-1:
        p[i] = r[i] + beta * p[i]

    rho0 = rho_new

输出结果 x

2. 并行主控流程伪代码

读入 CSR 矩阵 A 和向量 b
分配共享数组 x, r, p, Ap
分配线程数组与线程参数结构
初始化共享标量 alpha, beta, rho0, stop
初始化 barrier

for tid = 0 .. num_threads-1:
    填充 thread_data[tid]
    创建线程 worker(thread_data[tid])

for tid = 0 .. num_threads-1:
    等待线程结束

销毁 barrier
释放共享内存
输出结果 x

3. 并行工作线程伪代码

根据 tid 计算当前线程负责的数据范围

for i 属于本线程负责范围:
    x[i] = 0
    r[i] = b[i]
    p[i] = r[i]

计算本线程初始 rho0 局部和
执行一次全局规约
若初始残差已满足收敛条件:
    stop = 1

for iter = 0 .. max_iter-1 且 stop == 0:
    for i 属于本线程负责范围:
        Ap[i] = 稀疏矩阵第 i 行与 p 的乘积

    计算本线程 local_pAp
    执行一次全局规约
    由 0 号线程计算 alpha
    同步所有线程

    如果 stop == 1:
        跳出循环

    for i 属于本线程负责范围:
        x[i] = x[i] + alpha * p[i]
        r[i] = r[i] - alpha * Ap[i]

    计算本线程 local_rho_new
    执行一次全局规约
    由 0 号线程计算残差范数、beta，并判断是否收敛
    同步所有线程

    如果 stop == 1:
        跳出循环

    for i 属于本线程负责范围:
        p[i] = r[i] + beta * p[i]

    同步所有线程

4. 三种并行方案差异伪代码

block_reduce：

start, end = 连续分块划分结果

for i = start .. end-1:
    计算 Ap[i]

local_pAp = 本线程范围内 p[i] * Ap[i] 之和
partial_pAp[tid] = local_pAp

local_rho = 本线程范围内 r[i] * r[i] 之和
partial_rho[tid] = local_rho

由 0 号线程统一求和

block_lock：

start, end = 连续分块划分结果

for i = start .. end-1:
    计算 Ap[i]

local_pAp = 本线程范围内 p[i] * Ap[i] 之和
加锁
global_pAp += local_pAp
解锁

local_rho = 本线程范围内 r[i] * r[i] 之和
加锁
global_rho_new += local_rho
解锁

cyclic_reduce：

start = tid
end = n

for i = start; i < end; i += num_threads:
    计算 Ap[i]

local_pAp = 本线程负责下标上的 p[i] * Ap[i] 之和
partial_pAp[tid] = local_pAp

local_rho = 本线程负责下标上的 r[i] * r[i] 之和
partial_rho[tid] = local_rho

由 0 号线程统一求和

4.6 正确性验证思路

本实验的正确性验证主要从两个层面展开。

第一，验证算法结果是否正确。由于三个并行版本都严格沿用了串行版 CG 的更新公式，因此在同一输入数据下，理论上应与串行版得到一致或近似一致的结果。验证时可将串行程序 sparse 的输出作为参考，再分别比较三个并行程序的输出。

第二，验证并行实现是否正确。尽管并行版引入了线程划分与同步，但各线程只更新自己负责的数据区间，而全局点积量则通过规约或加锁处理，因此在设计上避免了逻辑数据竞争。实际验证时，应重点检查：

• 串行版与并行版输出结果是否一致
• 同一并行程序在不同线程数下输出是否一致
• 残差范数是否正常下降并满足收敛条件

按照这一思路，可以先在较小规模数据上完成正确性验证，再在正式实验规模下进行多次重复计时，从而保证后续性能分析建立在正确实现的基础上。

五、实现方法

5.1 串行程序实现

串行程序对应文件为 sparse.cpp。其实现主要包括矩阵与向量读入、spmv_csr 稀疏矩阵向量乘函数，以及 cg_serial 共轭梯度迭代函数。

其中，spmv_csr 负责按 CSR 存储格式完成矩阵向量乘；cg_serial 负责完成标准 CG 迭代流程。核心代码如下：

void spmv_csr (int n, const double* values, const int* row_ptr, const int* col_idx, const double* x, double* y)
{
    for (int i = 0; i < n; i++)
    {
        y[i] = 0.0;
    }
    for (int i = 0; i < n; i++)
    {
        for (int j = row_ptr[i]; j < row_ptr[i + 1]; j++)
        {
            y[i] += values[j] * x[col_idx[j]];
        }
    }
}

void cg_serial (int n, const double* values, const int* row_ptr, const int* col_idx, const double* b, double* x, int max_iter, double tol)
{
    double *r = (double*)malloc(n * sizeof(double));
    double *p = (double*)malloc(n * sizeof(double));
    double *Ap = (double*)malloc(n * sizeof(double));

    for (int i = 0; i < n; i++)
    {
        x[i] = 0.0;
        r[i] = b[i];
        p[i] = r[i];
    }

    double rho0 = 0.0;
    for (int i = 0; i < n; i++)
        rho0 += r[i] * r[i];

    for (int iter = 0; iter < max_iter; iter++)
    {
        spmv_csr(n, values, row_ptr, col_idx, p, Ap);

        double pAp = 0.0;
        for (int i = 0; i < n; i++)
            pAp += p[i] * Ap[i];

        double alpha = rho0 / pAp;
        for (int i = 0; i < n; i++)
        {
            x[i] += alpha * p[i];
            r[i] -= alpha * Ap[i];
        }

        double rho_new = 0.0;
        for (int i = 0; i < n; i++)
            rho_new += r[i] * r[i];

        if (sqrt(rho_new) < tol)
            break;

        double beta = rho_new / rho0;
        for (int i = 0; i < n; i++)
            p[i] = r[i] + beta * p[i];

        rho0 = rho_new;
    }
}

串行版的优点是结构清晰、逻辑完整、结果稳定，因此既可用于正确性基准，也可作为后续并行加速比的参考。

5.2 并行程序实现

三个并行程序分别对应：

• sparse_parallel_block_reduce.cpp
• sparse_parallel_block_lock.cpp
• sparse_parallel_cyclic_reduce.cpp

它们在整体结构上保持一致，主要由以下几部分组成：

• ThreadData 线程参数结构
• compute_range 任务划分函数
• cg_worker 工作线程函数
• cg_parallel 线程创建与回收函数

三个并行版本都使用多线程共享矩阵与向量数据，并通过 pthread_barrier_t 在各阶段之间同步。工作线程的主要执行过程为：

1. 初始化自己负责部分的 x、r、p
2. 计算局部初始残差平方和
3. 并行计算本地 Ap
4. 并行计算本地 p·Ap
5. 合并全局 p·Ap 后求 alpha
6. 并行更新 x 与 r
7. 并行计算本地 r·r
8. 合并全局 r·r 后求 beta
9. 并行更新 p

其中，三者的关键区别在于：

• block_reduce：连续分块 + 局部规约
• block_lock：连续分块 + 加锁规约
• cyclic_reduce：循环划分 + 局部规约

block_reduce 中连续分块的核心实现如下：

static void compute_range(int n, int num_threads, int tid, int& start, int& end)
{
    int base = n / num_threads;
    int extra = n % num_threads;

    start = tid * base + (tid < extra ? tid : extra);
    end = start + base + (tid < extra ? 1 : 0);
}

block_lock 中加锁规约的核心实现如下：

pthread_mutex_lock(data->mutex);
*(data->global_pAp) += local_pAp;
pthread_mutex_unlock(data->mutex);

pthread_mutex_lock(data->mutex);
*(data->global_rho_new) += local_rho_new;
pthread_mutex_unlock(data->mutex);

cyclic_reduce 中循环划分的核心实现如下：

static void compute_range(int n, int num_threads, int tid, int& start, int& end)
{
    start = tid;
    end = n;
}

for (int i = start; i < end; i += data->num_threads)
{
    data->Ap[i] = 0.0;
    for (int j = data->row_ptr[i]; j < data->row_ptr[i + 1]; j++)
    {
        data->Ap[i] += data->values[j] * data->p[data->col_idx[j]];
    }
}

从实现角度看，三个版本的数学求解流程相同，真正形成差异的是任务划分方式和全局点积规约方式。

5.3 数据生成程序

实验数据由 generate_input.cpp 生成。考虑到实验规模较大，生成器采用流式方式构造带状稀疏对称正定矩阵，并同步生成右端向量。

当前正式实验所用数据规模为：

• n = 5000000
• nnz = 34999988
• bandwidth = 3
• seed = 42

对应输入文件为：

• matrix.txt
• vector.txt

采用带状稀疏对称正定矩阵的原因是：

• 满足共轭梯度法对矩阵性质的要求
• 非零元个数与矩阵规模近似线性相关，便于构造大规模数据
• 能够控制文件大小和生成时间

5.4 PBS 脚本与运行方式

当前实验采用 PBS 提交作业，所有程序均在单节点共享内存环境下运行。实验目录中保留了以下脚本：

• 串行版本：sparse_1.pbs
• block_reduce：sparse_block_reduce_1/2/4/8/16.pbs
• block_lock：sparse_block_lock_1/2/4/8/16.pbs
• cyclic_reduce：sparse_cyclic_reduce_1/2/4/8/16.pbs

脚本中统一采用如下运行方式（以 sparse_block_reduce_8.pbs 为例）：

#!/bin/bash
#PBS -N sparse_br_8
#PBS -l nodes=1:ppn=8
#PBS -j oe

NAME="NiuTianhao"
ID="3024244288"
MATRIX_FILE="matrix.txt"
VECTOR_FILE="vector.txt"

cd $PBS_O_WORKDIR
procs=$(cat $PBS_NODEFILE | wc -l)

echo "=========================================="
echo "User: $NAME ($ID)"
echo "Program: sparse_parallel_block_reduce"
echo "Input: $MATRIX_FILE $VECTOR_FILE"
echo "Procs: $procs"
echo "Job: $PBS_JOBID"
echo "=========================================="

time ./sparse_parallel_block_reduce $MATRIX_FILE $VECTOR_FILE $procs >& run_block_reduce_8.log

其中：

• matrix.txt 为稀疏矩阵输入文件
• vector.txt 为右端向量输入文件
• $procs 为 PBS 分配的线程数

每组作业重复运行 10 次，再从输出文件中提取 real 时间汇总到 data.md。

六、结果分析

6.1 实验数据来源与统计口径

本节结果分析所用数据全部来源于服务器上重复提交 PBS 作业后的运行时间统计。时间取自 PBS 输出文件中的 real 字段，每组配置均重复运行 10 次，再整理到 data.md 中。

当前统计口径如下：

• 数据规模固定为 n = 5000000
• 串行基线为 sparse
• 并行方案包括 block_reduce、block_lock 和 cyclic_reduce
• 并行线程数取 1 / 2 / 4 / 8 / 16
• 每组实验均重复运行 10 次

需要说明的是，服务器环境下的 real 时间不仅反映程序计算时间，也会受到节点负载、文件缓存和作业调度状态的影响，因此不同批次之间存在一定波动。

6.2 正确性验证结果

根据程序结构和前期测试情况，串行版与三个并行版在相同输入下能够得到一致的求解结果，说明并行化没有改变 CG 的数学求解过程。

正确性验证主要基于以下事实：

• 三个并行版本都沿用了串行版的同一套更新公式
• 每个线程只更新自己负责的数据区间
• 全局点积量都经过规约或加锁处理

因此，可以认为当前实验中的四个版本在功能上是正确的，后续性能差异主要来自实现方式不同，而非算法错误。

6.3 平均运行时间统计

根据当前 data.md 中的数据，可得到各组平均运行时间如下：

方案	1线程/s	2线程/s	4线程/s	8线程/s	16线程/s
sparse	107.871	-	-	-	-
block_reduce	104.276	53.532	44.115	27.672	26.974
block_lock	109.282	58.365	39.575	28.974	26.788
cyclic_reduce	95.335	80.158	81.270	103.015	73.678

从表中可以看出：

• block_reduce 和 block_lock 随线程数增加整体明显变快
• block_reduce 在 8 线程后已接近性能上限
• block_lock 总体趋势正常，但 4 线程组有明显分层现象
• cyclic_reduce 基本没有表现出合理的并行加速趋势

6.4 加速比与并行效率分析

以串行基线 sparse 的平均时间 T_serial = 107.871s 为参考，可计算各并行方案的近似加速比和并行效率。

加速比与并行效率公式为：

$$S_p = \frac{T_{serial}}{T_p}$$ $$E_p = \frac{S_p}{p}$$

根据当前数据，可得到：

方案	2线程 S/E	4线程 S/E	8线程 S/E	16线程 S/E
block_reduce	2.015 / 1.008	2.445 / 0.611	3.898 / 0.487	3.999 / 0.250
block_lock	1.848 / 0.924	2.726 / 0.681	3.723 / 0.465	4.027 / 0.252
cyclic_reduce	1.346 / 0.673	1.327 / 0.332	1.047 / 0.131	1.464 / 0.092

这里需要说明两点：

1. 由于串行基线本身波动较大，且并行 1 线程版本与串行版实现细节不同，因此个别组会出现效率略大于 1 的现象。
2. 当前加速比和效率更适合用来比较趋势，而不应机械地按理想模型解释。

从趋势上看：

• block_reduce 与 block_lock 都获得了明显加速
• 8 线程是较明显的性能提升点
• 16 线程相对 8 线程提升已经较小，说明并行收益开始趋于饱和
• cyclic_reduce 的加速效果最差，几乎没有形成有效并行加速

进一步深度分析并行效率下降与加速饱和的核心原因：

1. 内存带宽瓶颈（访存受限 / Memory-Bound）： 稀疏矩阵向量乘（SpMV）涉及大量的间接内存访问（如 p[col_idx[j]]），其计算密度较低。这属于并行计算中典型的“内存受限”型应用。当线程数（如 16 线程）增加到一定规模时，多核并发高频的访存请求会触及甚至打满服务器节点内存总线的带宽物理上限，此时继续增加 CPU 核心已无法转化为有效的算力提升。
2. 迭代同步开销的累积： 共轭梯度法在单次迭代中需要进行多次全局规约求和与 pthread_barrier_wait 屏障同步。根据阿姆达尔定律（Amdahl’s Law），随着线程数暴增，线程间的长尾等待效应和锁同步周期会显著放大，这是导致规模较大时加速比趋于非线性的另一关键因素。

为更直观地展示各方案性能差异，以下给出各项指标的可视化对比图：

6.5 不同并行方案对比分析

1. block_reduce 与 block_lock

这两种方案都采用连续分块，因此主要差别来自规约方式。

从当前数据看：

• block_reduce 在 2 线程下略优于 block_lock
• block_lock 的 4 线程均值略快，但该组数据存在明显分层，稳定性较差
• 在 8 和 16 线程下，两者都收敛到约 27s 左右，性能非常接近

这说明在当前程序中，锁竞争虽然存在，但并没有成为压倒性的瓶颈。原因在于程序采用了合理的“粗粒度”锁定策略：block_lock 并没有在每次处理接矩阵元素求和时都去竞争全局锁，而是让每个线程先行局部独立完成 $O(N/T)$ 级别的大量点积运算任务，此后仅在最终规约环节进入短暂的互斥临界区。既然进锁的频次极低且执行迅速，这部分通信开销就被庞大的本地计算量所彻底摊派抵消了。

2. block_reduce 与 cyclic_reduce

这组对比最具有代表性，因为两者都采用局部规约，主要差异只来自任务划分方式。

实验结果表明：

• block_reduce 能随着线程数增加稳定加速
• cyclic_reduce 即使在 16 线程下仍明显慢于 block_reduce
• cyclic_reduce_8 甚至比 1 线程还慢

这说明在当前问题中，连续分块不仅具有压倒性优势，而且从并行设计的机理上看必然优越于循环划分法。深度原因主要包括：

• 无谓的负载均衡开销：结合 5.3 节代码可知，实验采用的是带状稀疏矩阵（非零元聚集且分布规整，各行计算量基本一致）。对于这种计算负载天然均衡的任务，连续分块已经能够实现很好的负载均衡。循环划分（Cyclic）通常用于缓解非均匀任务带来的负载失衡，但在本实验中不仅未能带来明显收益，反而引入了额外的开销。
• 破坏缓存局部性（Cache Locality）：连续分块能够保证线程访问连续的内存空间，从而充分利用 CPU 缓存，保持较高的缓存命中率。相比之下，循环划分跳跃式地访问数组和矩阵元素，破坏了空间局部性，导致缓存未命中（Cache Miss）率显著上升，增加了访存延迟。
• 严重的伪共享（False Sharing）问题：在多线程场景下，cyclic_reduce 暴露出底层的性能瓶颈。由于相邻的 double 变量极可能位于同一个 CPU 缓存行（Cache Line）中，多线程在交错更新 x、r、p、Ap 等数组时，会频繁修改同一缓存行内的不同数据。这触发了缓存一致性协议，导致缓存行在不同核心间不断失效和同步（缓存乒乓效应），带来了巨大的额外开销，最终导致性能严重退化。

3. 关于异常波动的分析

当前数据中存在多组明显分层现象，例如：

• 串行 sparse 前后两档明显不同
• block_lock_4 前 6 次和后 4 次形成两档
• cyclic_reduce_8 前 8 次和后 2 次形成两档

这说明实验环境本身存在一定波动，可能原因包括：

• 节点负载变化
• 文件缓存状态变化
• PBS 调度和 I/O 开销波动

因此，在分析结果时不应只看某一次时间，而应结合 10 次均值与整体波动情况进行判断。

4. 当前阶段的总体结论

结合已有结果，可以得到以下结论：

1. 连续分块方案整体优于循环划分方案。
2. 在当前实现中，局部规约与加锁规约的性能差距并不大。
3. block_reduce 是最稳定、最符合预期的并行方案。
4. cyclic_reduce 的性能明显异常，主要原因很可能是较差的局部性和较严重的伪共享，再叠加环境波动影响。
5. 对于当前规模 n = 5000000，程序在 8 线程以后已接近平台上的共享内存并行上限，继续增加线程带来的收益有限。

七、个人总结

这次实验不仅是对 CSR 稀疏矩阵和共轭梯度法（CG）代码的堆彻，更像是一次对计算机底层硬件特性的摸底。动手写完三种 Pthread 并行方案并出完数据后，我最大的感触是：并行计算真不是简单地“多加几个线程程序就一定能变快”。想要有好的加速效果，关键得看任务怎么分、锁怎么加，以及有没有照顾到 CPU 缓存的物理特性。

从实验数据能很直观地看到，原本以为能把任务分得更均匀的循环划分（Cyclic），反倒因为破坏了缓存局部性、引发了严重的缓存乒乓效应（伪共享），最终在多核下被最基础的连续分块（Block）全面碾压。这说明顺着内存空间局部性的特点写代码，比生硬地追求“绝对平均分工”更重要。同时，block_lock 和 block_reduce 性能差不多，也让我意识到只要把加锁的粒度控制好（比如只在最后归约时加极少次数的锁），锁带来的开销其实能被庞大的计算量轻松摊平。

此外，由于集群节点会有其他人跑任务，我测数据时也切实感觉到了性能波动，认识到跑实验得多测几次看平均趋势才靠谱。值得一提的是，在这次实验的数据图表生成以及底层硬件理论（如伪共享）的梳理排查中，我主动借助了大模型工具作为辅助，这有效提升了我的科研数据处理效率。总之，这次实验让我把课本上的多线程编程、锁机制、高速缓存等知识真正和代码串联了起来，收获非常扎实。

附录

A.1 源码与可执行文件

• sparse.cpp
• sparse_parallel_block_reduce.cpp
• sparse_parallel_block_lock.cpp
• sparse_parallel_cyclic_reduce.cpp
• generate_input.cpp

A.2 数据文件与统计文件

• matrix.txt
• vector.txt
• data.md
• 流程图.png

A.3 PBS 脚本文件

• sparse_1.pbs
• sparse_block_reduce_1.pbs
• sparse_block_reduce_2.pbs
• sparse_block_reduce_4.pbs
• sparse_block_reduce_8.pbs
• sparse_block_reduce_16.pbs
• sparse_block_lock_1.pbs
• sparse_block_lock_2.pbs
• sparse_block_lock_4.pbs
• sparse_block_lock_8.pbs
• sparse_block_lock_16.pbs
• sparse_cyclic_reduce_1.pbs
• sparse_cyclic_reduce_2.pbs
• sparse_cyclic_reduce_4.pbs
• sparse_cyclic_reduce_8.pbs
• sparse_cyclic_reduce_16.pbs