并行计算与性能优化实践

1、对于台式机、笔记本电脑或手机，其系统的理论并行处理能力与串行处理能力相比如何？其中存在哪些类型的并行硬件？

很难穿透营销噱头找到真实规格。大多数设备（包括手持设备）都有多核处理器和至少一个集成图形处理器。除了非常老旧的硬件，台式机和笔记本电脑都有一些向量计算能力。

2、你有一个图像处理应用程序，每天需要处理1000张图像，每张图像大小为4兆字节（MiB，2的20次方或1048576字节）。串行处理每张图像需要10分钟。你的集群由多核节点组成，每个节点有16个核心，总主内存存储为16吉字节（GiB，2的30次方字节，即1024兆字节）。哪种并行处理设计最适合处理这个工作负载？

单计算节点上的线程化以及向量化。处理时每次处理16张图像所需内存为 4 MiB × 16 = 64 MiB，远低于集群中每个节点的 16 GiB 内存。串行处理总时间为 10分钟 × 1000 = 10000 分钟，换算为小时约为 167 小时（10000 ÷ 60 ≈ 167）。在 16 个核心上并行处理时间约为 10000 ÷ 16 = 625 分钟，约 10.4 小时（625 ÷ 60 ≈ 10.4）。向量化可将时间缩短至 5 小时以内。

3、你有一个图像处理应用程序，每天需要处理1000张图像，每张图像大小为4兆字节（MiB，2的20次方或1048576字节）。串行处理每张图像需要10分钟。你的集群由多核节点组成，每个节点有16个核心，主内存存储总量为16吉字节（GiB，2的30次方字节，即1024兆字节）。现在客户需求增加了10倍。你的设计能处理这种情况吗？你需要做出哪些改变？

原设计在客户需求增加10倍后可能仍可行，但处理时间会变为100分钟。此时可能需要考虑采用 消息传递 或 分布式计算 。

4、你有一个在研究生阶段开发的波高模拟应用程序。它是一个串行应用程序，并且由于它原本仅被计划作为你博士论文的基础，你没有采用任何软件工程技术。现在你计划将其作为一个可供许多研究人员使用的工具的起点。你的团队还有另外三名开发人员。你会在这个项目计划中包含哪些内容？

项目计划内容

项目计划可包含以下内容：

设置版本控制以跟踪代码变化开发测试套件来确保应用程序结果的准确性清理现有代码，提升代码质量和可移植性确定可用计算资源的能力、应用程序的需求和性能要求，可通过系统基准测试和性能分析来实现根据内核分析结果，规划并行化例程的任务并实施更改实施更改后，将增量更改提交到版本控制系统建立团队开发流程，进行有效的项目管理，包括任务管理和范围控制

5、使用CTest创建一个测试


1. 前提条件：需要安装MPI、CMake和ndiff。MPI使用OpenMPI 4.0.0，CMake使用3.13.3（包含CTest），GCC编译器使用版本8。在Mac上使用Homebrew，在Ubuntu Linux上使用Apt和Synaptic进行安装，确保从libopenmpi-dev获取开发头文件。ndiff需从[https://www.math.utah.edu/~beebe/software/ndiff/](https://www.math.utah.edu/~beebe/software/ndiff/)下载工具，运行`./configure`、`make`和`make install`手动安装。

2. 创建源文件：制作两个源文件，如`TimeIt.c`和`MPITimeIt.c`。`TimeIt.c`是C程序，`MPITimeIt.c`是MPI程序。

3. 在CMakeLists.txt中添加：
   ```cmake
   enable_testing()
   add_test(NAME make WORKING_DIRECTORY ${CMAKE_BINARY_DIRECTORY} COMMAND ${CMAKE_CURRENT_SOURCE_DIR}/build.ctest)
   ```

4. 添加build.ctest文件：
   ```sh
   #!/bin/sh
   if [ -x $0 ]
   then
       echo "PASSED - is executable"
   else
       echo "Failed - ctest script is not executable"
       exit -1
   fi
   ```

5. 运行测试：
   ```bash
   mkdir build && cd build
   cmake ..
   make
   make test
   ```
   或者：
   ```bash
   ctest
   ```
   也可以使用以下命令获取失败测试的输出：
   ```bash
   ctest --output-on-failure
   ```

6、在你选择的一个小型应用程序上运行Valgrind

可以按以下步骤在小型应用程序上运行 Valgrind：

准备运行环境，可在不同计算机、旧版 macOS 上运行，或创建虚拟机或 Docker 镜像。若使用 GCC 编译器，其与 Valgrind 配合效果好；若使用英特尔编译器，编译时不进行向量化以避免向量指令警告。运行程序时，将 valgrind 命令置于可执行程序名前。对于 MPI 作业， valgrind 命令放在 mpirun 之后、可执行程序名之前。

例如，对于示例代码，先使用 gcc -g -o test test.c 编译代码，然后用 valgrind --leak-check=full ./test 2 运行。

7、计算一个系统的理论性能。在计算中包括峰值浮点运算次数、内存带宽和机器平衡。假设已知该系统的峰值浮点运算次数为 236.8 GFlops/s，内存带宽为 34.1 GiB/s。

理论机器平衡（MBT）的计算如下：

因为 1 字节等于 8 比特，所以将内存带宽单位换算为字节后，理论机器平衡 MBT = 峰值浮点运算次数（FT） / 内存带宽（BT） × (8 字节/字)

即：


MBT = 236.8 GFlops/s / 34.1 GiB/s × (8 字节/字) = 56 Flops/字

该系统的峰值浮点运算次数为 236.8 GFlops/s ，内存带宽为 34.1 GiB/s ，理论机器平衡为 56 Flops/字 。

8、从https://bitbucket.org/berkeleylab/cs – roofline – toolkit.git下载Roofline工具包，并测量你所选系统的实际性能。

可通过以下步骤实现：

先使用命令
git clone https://bitbucket.org/berkeleylab/cs-roofline-toolkit.git
下载Roofline工具包；

若要测量NVIDIA V100 GPU的性能：
– 输入
cd cs-roofline-toolkit/Empirical_Roofline_Tool-1.1.0
– 再执行
cp Config/config.voltar.uoregon.edu Config/config.V100_gpu
– 编辑 Config/config.V100_gpu 中的设置（如：
ERT_RESULTS Results.V100_gpu ERT_PRECISION FP64 ERT_NUM_EXPERIMENTS 5
等）；
– 运行
tests ./ert Config/config.V100_gpu
命令；
– 最后查看
Results.config.V100_gpu/Run.001/roofline.ps
文件；

若要测量AMD Vega20 GPU的性能：
– 执行
cp Config/config.odinson-ocl-fp64.01 Config/config.Vega20_gpu
– 编辑 Config/config.Vega20_gpu 中的设置（如：
ERT_RESULTS Results.Vega20_gpu ERT_CFLAGS -O3 -x c++ -std=c++11 -Wno-deprecated-declarations -I<OpenCL头文件路径> ERT_LDLIBS -L<OpenCL库文件路径> -lOpenCL
等）；
– 运行
tests ./ert Config/config.Vega20_gpu
– 查看
Results.config.Vega20_gpu/Run.001/roofline.ps
文件中的输出。

9、确定一个小型应用程序的平均处理器频率和能耗


可以使用 `likwid` 工具套件中的 `likwid-powermeter` 命令行工具查看处理器频率和功率统计信息，`likwid-perfctr` 工具也会在摘要报告中报告部分统计信息；还可以使用 Intel® Power Gadget，它有适用于 Mac 和 Windows 的版本，Linux 版本功能更有限，能以图表形式展示频率、功率、温度和利用率。CLAMR 迷你应用正在开发的 PowerStats 库可在应用程序内跟踪能量和频率，并在运行结束时报告，目前该库在 Mac 上可使用 Intel Power Gadget 库接口工作，Linux 系统的类似功能也在开发中，应用代码只需添加 `powerstats_init(); powerstats_sample(); powerstats_finalize();` 这几个调用。

10、确定一个应用程序使用了多少内存

首先，通过 top 或 ps 命令获取进程 ID。然后使用以下命令之一来跟踪内存使用情况：

watch -n 1 "grep VmRSS /proc/<pid>/status" watch -n 1 "ps <pid>" top -s 1 -p <pid>

此外，CLAMR 中的 MemSTATS 库提供了四个不同的内存跟踪调用，可将其集成到程序中，以查看不同阶段的内存使用情况。

11、编写一个C语言的二维内存分配器，使其内存布局与Fortran相同。

下面是给定的【文本内容】：

假设在 Fortran 中按 x(j,i) 访问数组，在 C 中按 x[i][j] 访问。创建宏


#define x(j,i) x[i - 1][j - 1]

代码就能使用 Fortran 数组表示法。可将二维内存分配器中 i 和 j 、 imax 和 jmax 互换。代码如下：


#include <stdlib.h>
#include "malloc2Dfort.h"

double **malloc2Dfort(int jmax, int imax) {
    double **x = (double **)malloc(imax*sizeof(double *) + imax*jmax*sizeof(double));
    x[0] = (double *)(x + imax);
    for (int i = 1; i < imax; i++) {
        x[i] = x[i - 1] + jmax;
    }
    return(x);
}

12、AVX – 512向量单元会如何改变流三元组的ECM模型？

使用AVX-256向量单元的ECM模型性能分析

在使用AVX-256向量单元的情况下，可在 1个周期内处理所有所需的浮点运算 。当使用AVX-512时，虽然仍需 1个周期 ，但 只有一半的向量单元处于忙碌状态 。若存在可完成两倍工作量的任务，AVX-512具备相应的执行能力。由于计算操作时间 TOL仍为1个周期 ，因此 性能不会发生任何变化 。

13、火山灰羽流中的云碰撞模型针对距离在1毫米以内的粒子启动。请编写空间哈希实现的伪代码。此操作的复杂度阶是多少？

伪代码如下：


for all particles, ip, in NParticles {
    for all particles, jp, in Adjacent_Buckets {
        if (distance between particles < 1mm) {
            perform collision or interaction calculation
        }
    }
}

此操作的复杂度阶是Θ(N)。

14、一个波浪模拟代码使用自适应网格细化（AMR）网格来更好地细化海岸线。模拟要求是记录指定位置（浮标和海岸设施所在位置）的波浪高度随时间的变化。由于单元格在不断细化，如何实现这一点？

创建完美空间哈希

创建一个完美的空间哈希，其箱大小与最小单元格相同。将单元格索引存储在单元格下方的箱中。为每个站点计算箱，并从箱中获取单元格索引。

15、对第4.3节多材料代码（https://github.com/LANL/MultiMatTest.git）中的循环进行自动向量化实验。添加向量化和循环报告标志，查看编译器给出的信息。

可按照要求对指定代码中的循环进行自动向量化实验，添加相应标志后，编译器会给出关于循环是否向量化、向量依赖情况、估计加速比等信息。

如使用 Intel 编译器，可能会显示以下信息：

循环是否被向量化是否存在向量依赖向量长度估计加速比

使用 GCC 9 及以上版本编译器也能对部分代码进行向量化。

16、为一个归约操作编写一个高级 OpenMP 版本

在高级 OpenMP 中，手动划分数据。数据分解在代码的第 6 – 9 行完成。线程 0 在第 13 行分配 xmax_thread 共享数据数组。第 18 – 22 行找出每个线程的最大值，并将结果存储在 xmax_thread 数组中。然后，在第 26 – 30 行，一个线程找出所有线程中的最大值。代码示例如下：


#include <stdlib.h>
#include <float.h>
#include <omp.h>
double array_max(double* restrict var, int ncells)
{
    int nthreads = omp_get_num_threads();
    int thread_id = omp_get_thread_num();
    int tbegin = ncells * ( thread_id ) / nthreads;
    int tend = ncells * ( thread_id + 1 ) / nthreads;
    static double xmax;

17、为什么在进行数据发送和接收时，不能仅在接收时进行阻塞操作，而在使用打包或数组缓冲区方法进行数据发送时采用非阻塞操作呢？

使用打包或数组缓冲区的版本会调度发送，但在数据复制或发送完成之前就返回了。

MPI_Isend标准规定：

“在非阻塞发送操作调用后，直到发送完成，发送方不应修改发送缓冲区的任何部分”。

打包和数组版本在通信后会释放缓冲区，所以这些版本可能会在缓冲区数据被复制之前就删除缓冲区，导致程序崩溃。

为了安全起见，在删除缓冲区之前必须检查发送状态。

18、在幽灵交换的向量类型版本中，进行阻塞接收是否安全？如果只进行阻塞接收，有什么优点？

向量版本直接从原始数组发送数据，而非复制数据，比分配缓冲区并释放的版本更安全。若只进行阻塞接收，通信速度会更快。

19、尝试在其中一个幽灵交换例程中，用 MPI_ANY_TAG 替换显式标签。这样做是否可行？是否会更快？使用显式标签有什么优势？

使用 MPI_ANY_TAG 作为标签参数可行，可能会稍快，但提升幅度可能小到难以测量。使用显式标签可额外检查是否接收到正确的消息。

20、在可以访问电源硬件计数器的系统上，使用likwid性能工具获取CloverLeaf应用程序的CPU功率需求。

可按以下步骤操作：

从包管理器安装 likwid 或使用以下命令安装：

bash git clone https://github.com/RRZE-HPC/likwid.git cd likwid edit config.mk make make install

使用 sudo modprobe msr 启用 MSR。

运行以下命令：

bash likwid-perfctr -C 0-87 -g MEM_DP ./clover_leaf

运行该命令后，在输出结果中可以找到 Power [W] STAT 对应的数值，此数值即为 CPU 功率需求。

例如运行结果中 Power [W] STAT 的值为 279.9804 （平均为 3.1816 ）。

21、有一个图像分类应用程序，将每个文件传输到 GPU 需要 5 毫秒，处理需要 5 毫秒，再传输回来需要 5 毫秒。在 CPU 上，每张图像的处理时间为 100 毫秒。有一百万张图像需要处理，CPU 有 16 个处理核心。GPU 系统处理这些工作会更快吗？

不会，CPU 处理时间为 6250 秒，GPU 处理时间为 15000 秒，GPU 系统处理时间约是 CPU 的 2.5 倍。

22、问题中GPU的传输时间基于第三代PCI总线。如果能使用第四代PCI总线，设计会有怎样的变化？使用第五代PCI总线呢？对于图像分类，不需要返回修改后的图像，这会如何改变计算结果？

第四代PCI总线速度是第三代的两倍，所需时间为：

(2.5 ms+5 ms+2.5 ms)×1,000,0001,000=10,000 s(2.5 ms+5 ms+2.5 ms)×1,000,0001,000=10,000 s

第五代PCI总线速度是原第三代的四倍，所需时间为：

(1.25 ms+5 ms+1.25 ms)×1,000,0001,000=7,500 s(1.25 ms+5 ms+1.25 ms)×1,000,0001,000=7,500 s

若不需要返回修改后的图像，计算时可去掉返回时间。

23、对于你的独立 GPU（如果没有，则为 NVIDIA GeForce GTX 1060），你可以运行多大规模的 3D 应用程序？假设每个单元格有 4 个双精度变量，并且使用一半的 GPU 内存，以便为临时数组留出空间。如果使用单精度，情况会如何变化？

对于双精度，可运行 465×465×465 的 3D 网格应用程序；对于单精度，可运行 586×586×586 的 3D 网格应用程序，使用单精度时分辨率提高了 25% 。

24、在本地 GPU 开发系统上运行 OpenACC/StreamTriad 和/或 OpenMP/StreamTriad 目录中的流三元组示例。你可以在 https://github.com/EssentialsofParallelComputing/Chapter11 找到这些目录。

要完成该任务，首先需确保本地 GPU 开发系统环境已正确配置，具备运行相关代码的条件。然后通过命令行工具或版本控制工具，从指定链接（ https://github.com/EssentialsofParallelComputing/Chapter11 ）克隆代码仓库到本地。

接着，进入 OpenACC/StreamTriad 或 OpenMP/StreamTriad 目录，根据不同编译器的编译指令和依赖库要求，编写合适的编译脚本进行编译。

编译成功后，执行生成的可执行文件，即可运行流三元组示例。

25、修改 OpenMP 数据区域映射，以反映内核中数组的实际使用情况。

将数据分配和释放部分，分别修改为使用 #pragma omp target enter data 和 #pragma omp target exit data 指令。具体修改如下：

原分配部分 13 #pragma omp target 及后续 malloc 操作，修改为
c 13 #pragma omp target enter data map(alloc:a[0:nsize], b[0:nsize], c[0:nsize])

原释放部分 42 #pragma omp target 及后续 free 操作，修改为
c 36 #pragma omp target exit data map(delete:a[0:nsize], b[0:nsize], c[0:nsize])

完整修改后的代码在 Stream_par7.c 文件中。

26、用OpenMP实现质量求和示例，假设有 `ncells` 个单元格，每个单元格有类型 `celltype` 、高度 `H` 、 `x` 方向的长度 `dx` 和 `y` 方向的长度 `dy` ，当单元格类型为 `REAL_CELL` 时，将该单元格的质量（高度乘以 `x` 方向长度乘以 `y` 方向长度）累加到总和中，编写代码实现该功能。


#include <stdio.h>
#define REAL_CELL 1

double mass_sum(int ncells, int* restrict celltype, double* restrict H, double* restrict dx, double* restrict dy) {
    double summer = 0.0;
    #pragma omp target teams distribute parallel for simd reduction(+:summer)
    for (int ic = 0; ic < ncells; ic++) {
        if (celltype[ic] == REAL_CELL) {
            summer += H[ic] * dx[ic] * dy[ic];
        }
    }
    return summer;
}

27、在 CUDA 流三元组示例中，将主机内存分配改为使用固定内存，性能有提升吗？

要使用固定内存，在主机端内存分配时，将 malloc 替换为 cudaHostMalloc ，释放内存时将 free 替换为 cudaFreeHost 。与不使用固定内存的情况进行性能比较，使用固定内存时，数据传输时间至少快两倍。

28、在GPU设备上使用SYCL初始化数组a和b

以下代码展示了在GPU上初始化数组 a 和 b 的版本：


#include <iostream>
#include <vector>
#include <chrono>
#include <CL/sycl.hpp>

namespace sycl = cl::sycl;

int main() {
    const size_t nsize = 1000; // 假设数组大小为1000
    std::vector<double> a(nsize);
    std::vector<double> b(nsize);
    std::vector<double> c(nsize);

    auto t1 = std::chrono::high_resolution_clock::now();

    // 选择GPU设备
    sycl::queue Queue(sycl::gpu_selector{}); 
    const double scalar = 3.0;

    sycl::buffer<double, 1> dev_a { a.data(), sycl::range<1>(a.size()) };
    sycl::buffer<double, 1> dev_b { b.data(), sycl::range<1>(b.size()) };
    sycl::buffer<double, 1> dev_c { c.data(), sycl::range<1>(c.size()) };

    Queue.submit([&](sycl::handler& CommandGroup) {
        auto a_acc = dev_a.get_access<sycl::access::mode::write>(CommandGroup);
        auto b_acc = dev_b.get_access<sycl::access::mode::write>(CommandGroup);
        auto c_acc = dev_c.get_access<sycl::access::mode::write>(CommandGroup);

        CommandGroup.parallel_for<class StreamTriad>( sycl::range<1>{nsize}, [=] (sycl::id<1> it) {
            a_acc[it] = 1.0;
            b_acc[it] = 2.0;
            c_acc[it] = -1.0;
        });
    });

    Queue.wait();

    Queue.submit([&](sycl::handler& CommandGroup) {
        auto a_acc = dev_a.get_access<sycl::access::mode::read>(CommandGroup);
        auto b_acc = dev_b.get_access<sycl::access::mode::read>(CommandGroup);
        auto c_acc = dev_c.get_access<sycl::access::mode::write>(CommandGroup);

        CommandGroup.parallel_for<class StreamTriad>( sycl::range<1>{nsize}, [=] (sycl::id<1> it) {
            c_acc[it] = a_acc[it] + scalar * b_acc[it];
        });
    });

    Queue.wait();

    auto t2 = std::chrono::high_resolution_clock::now();
    double time1 = std::chrono::duration_cast<std::chrono::duration<double> >(t2 - t1).count();
    std::cout << "Runtime is " << time1*1000.0 << " msecs " << std::endl;

    return 0;
}

注意：

原答案中使用了 sycl::cpu_selector{} ，这里修改为 sycl::gpu_selector{} 以确保使用GPU设备。避免了变量名重复，将 a 、 b 、 c 访问器重命名为 a_acc 、 b_acc 、 c_acc 。补充了必要的头文件和 main 函数结构，使代码可以独立编译运行。