CS149ParallelComputing_Assignments

repo链接

照着23年的Assignments写的

  • Assignment 1: Analyzing Parallel Program Performance on a Quad-Core CPU
  • Assignment 2: Scheduling Task Graphs on a Multi-Core CPU
  • Assignment 3: A Simple Renderer in CUDA
  • Assignment4:NanoGPT149

个人体感是第一个实验虽然是最简单的,但我从第一个实验中学到最多

Assignment1

环境:

  • CPU: 11th Gen Intel(R) Core(TM) i5-11400H @ 2.70GHz
  • 操作系统:wls2(Ubuntu 24.04.1 LTS) 限制使用 8个 (virtual)cores, 12GB 内存. 参考 how to modify cpu cores in wls2

该PA一共要求要完成6个小实验

Prog1 Parallel Fractal Generation Using Threads

生成两张图,右图是左经过放大后得来的,其中较白像素的计算量大于较黑像素的计算量
image-20260624230657266

试验具体要求是让我们使用多线程的方式,加速图片的生成,并检查加速比(多线程版本相对于单线程版本)是否随着线程数的增加而线性地增加

第一个实现方式,均等地将图片划分成几个块,分给不同的线程工作

void workerThreadStart(WorkerArgs * const args) {
    double startTime = CycleTimer::currentSeconds();
    // 1. method 1 , cannot achieve high efficiency because of loading imbalance
    int rowsPerThread = args->height / args->numThreads;
    int startRow = args->threadId * rowsPerThread;
    int totalRows = rowsPerThread;
    // avoid remaining rows not calculated
    if (args->threadId + 1 == args->numThreads) {
        totalRows = args->height - startRow;
    }
    mandelbrotSerial(args->x0, args->y0, args->x1, args->y1, 
    args->width, args->height, startRow, totalRows, args->maxIterations, args->output);
    double endTime = CycleTimer::currentSeconds();
    printf("Thread %d completed in %.3f ms\n", 
           args->threadId, (endTime - startTime) * 1000);
}

结果如下,随着线程数量增多,加速比并没有线性地增加,拿VIew1举例来说,当我们的thread增加到8时,达到加速比 为4.04 ,远小于预期的8:

VIEW 1:
Threads | Serial Time | Thread Time | Speedup
--------|-------------|-------------|--------
   1 |      443.88 |      441.42 |    1.01
   2 |      429.90 |      217.44 |    1.98
   3 |      431.69 |      273.03 |    1.58
   4 |      447.57 |      184.53 |    2.43
   5 |      470.56 |      188.15 |    2.50
   6 |      441.93 |      136.66 |    3.23
   7 |      445.93 |      131.33 |    3.40
   8 |      444.61 |      110.16 |    4.04
  16 |      444.29 |       77.78 |    5.71

VIEW 2:
Threads | Serial Time | Thread Time | Speedup
--------|-------------|-------------|--------
   1 |      273.06 |      265.42 |    1.03
   2 |      263.66 |      155.34 |    1.70
   3 |      275.75 |      124.29 |    2.22
   4 |      280.07 |      110.14 |    2.54
   5 |      273.51 |       90.61 |    3.02
   6 |      274.46 |       79.32 |    3.46
   7 |      271.91 |       71.36 |    3.81
   8 |      273.92 |       63.18 |    4.34
  16 |      273.22 |       40.52 |    6.74

为什么有这样的现象?因为图片在纵向上的负载分配极不均衡,拖累了整体的加速比,为了证明这个观点,可以把每一个线程的执行时间打印出来。如下所示,对于VIew1来说“中间”的线程执行时间最长,而“上下两端”的线程执行时长很短

8 threads for view1
Thread 0 completed in 7.596 ms    // 非常短
Thread 7 completed in 7.823 ms      // 非常短
Thread 1 completed in 38.047 ms
Thread 6 completed in 38.263 ms
Thread 2 completed in 74.456 ms
Thread 5 completed in 76.675 ms
Thread 3 completed in 110.445 ms  // 非常长
Thread 4 completed in 114.860 ms  // 非常长

8 threads for view2   // view2 的负载比view1好些,但也存在肉眼可见的不均衡
Thread 5 completed in 27.754 ms
Thread 2 completed in 28.456 ms
Thread 6 completed in 32.051 ms
Thread 4 completed in 32.241 ms
Thread 7 completed in 35.941 ms
Thread 1 completed in 40.988 ms
Thread 3 completed in 45.732 ms
Thread 0 completed in 61.845 ms   // 最顶端的线程负载最重,符合直观预期

修改分配方式,实现2:对图像的每一行交叉相错地分配给某个线程,这样每个线程的负载能达到更好的平衡

void workerThreadStart(WorkerArgs * const args) {
    double startTime = CycleTimer::currentSeconds();
    // 2. method2 Interleaved row assignment for better load balancing
    // Thread i processes rows i, i+numThreads, i+2*numThreads, etc.
    for (int row = args->threadId; row < args->height; row += args->numThreads) {
        mandelbrotSerial(args->x0, args->y0,
                        args->x1, args->y1,
                        args->width, args->height, row, 
                        1, args->maxIterations, args->output);
    }
    double endTime = CycleTimer::currentSeconds();
    printf("Thread %d completed in %.3f ms\n", 
           args->threadId, (endTime - startTime) * 1000);
}

结果如下所示,可以看到当线程数为8时,我们对VIEW1的加速比达到7.31,非常接近理论值8。而且有意思的是,当我们把线程数增加到16时,加速比下降了!对比实现1,8 -> 16的线程数增加会增加加速比。为什么有这样的区别?我个人认为实现二任务划分的颗粒度很小,线程越多则导致系统用于通信、同步、缓存同步的开销就越多,导致了加速比下降。但对于实现一,我们的颗粒度较大,如果增加线程,其负载均衡的正收益大于通信、同步带来的负收益。(Lecture5 的理论课程应证了实验现象)

VIEW 1 :
Threads | Serial Time | Thread Time | Speedup
--------|-------------|-------------|--------
   1 |      441.05 |      441.50 |    1.00
   2 |      434.64 |      221.71 |    1.96
   3 |      442.01 |      150.47 |    2.94
   4 |      443.59 |      114.21 |    3.88
   5 |      442.21 |       91.33 |    4.84
   6 |      442.64 |       78.35 |    5.65
   7 |      435.25 |       67.65 |    6.43
   8 |      441.04 |       60.35 |    7.31
  16 |      441.06 |       62.60 |    7.05

VIEW 2 :
Threads | Serial Time | Thread Time | Speedup
--------|-------------|-------------|--------
   1 |      268.30 |      273.04 |    0.98
   2 |      271.73 |      138.59 |    1.96
   3 |      272.25 |       93.52 |    2.91
   4 |      271.67 |       70.59 |    3.85
   5 |      273.56 |       58.18 |    4.70
   6 |      273.52 |       49.04 |    5.58
   7 |      273.04 |       42.22 |    6.47
   8 |      265.13 |       39.06 |    6.79
  16 |      267.31 |       40.03 |    6.68

再打印他们的各线程的执行时间,明显地看到其负载相对于实现1来说更加均衡了:

view1:
Thread 5 completed in 56.712 ms
Thread 6 completed in 57.062 ms
Thread 2 completed in 58.365 ms
Thread 4 completed in 58.699 ms
Thread 3 completed in 58.780 ms
Thread 1 completed in 59.752 ms
Thread 7 completed in 82.346 ms
Thread 0 completed in 85.910 ms
view2:
Thread 2 completed in 37.255 ms
Thread 6 completed in 34.747 ms
Thread 5 completed in 35.025 ms
Thread 1 completed in 35.274 ms
Thread 4 completed in 35.364 ms
Thread 3 completed in 35.932 ms
Thread 2 completed in 36.120 ms
Thread 7 completed in 35.609 ms
Thread 0 completed in 38.798 ms

但无论如何都做不到随着线程数增加加速比“线性地”增加,因为多线程必然带来通信、同步、缓存一致性等开销。

Prog2 Vectorizing Code Using SIMD Intrinsics

虽说标题是“Intrinsics”,但实际上还是在软件层模拟了SIMD的运行方式,其目的是让我们能“可视化”地体会真实SIMD的实现。

具体代码实现就不贴了,可以参考git仓库,有两个注意点:

  1. 当迭代次数不是 vector_width的整数倍数时,需要做特殊处理
  2. 查看SIMD函数的实现方式,几乎都有一个mask,表示某个通道是否要进行计算,所有通道的判断结果越趋同则表示通道利用率越高,同样表示divergent execution的程度越低,这在真实的SIMD场景中(ISPC以及CUDA)往往取得更好的性能

试验结果如下,vector_width越大,通道利用率越低,这是符合预期的

VECTOR_WIDTH | Vector Utilization 
-------------|-------------
   2         |    87.6%
   4         |    82.3%
   8         |    79.6%
   16        |    78.3%

Prog3 Parallel Fractal Generation Using ISPC

part1

试验结果:

view1:
[mandelbrot serial]:            [217.952] ms
Wrote image file mandelbrot-serial.ppm
[mandelbrot ispc]:              [44.823] ms
Wrote image file mandelbrot-ispc.ppm
                                (4.86x speedup from ISPC)

view2 :
[mandelbrot serial]:            [134.371] ms
Wrote image file mandelbrot-serial.ppm
[mandelbrot ispc]:              [32.859] ms
Wrote image file mandelbrot-ispc.ppm
                                (4.09x speedup from ISPC)

理论上,我可以取得最大8的加速比(因为ISPC编译器生成8-wide AVX2的向量指令),但实际上两个图都只取得了4点多的加速比。

为什么?

part1的ispc代码只会利用一个核,所以它的加速比只来自于底层的AVX2指令(SIMD)。下面的ispc代码,从上到下地遍历图像,且会按照8像素一组的方式进行计算,每一组都同时计算这8个像素的值(调用mandel)。

export void mandelbrot_ispc(uniform float x0, uniform float y0, 
                            uniform float x1, uniform float y1,
                            uniform int width, uniform int height, 
                            uniform int maxIterations,
                            uniform int output[])
{
    float dx = (x1 - x0) / width;
    float dy = (y1 - y0) / height;

    foreach (j = 0 ... height, i = 0 ... width) {
            float x = x0 + i * dx;
            float y = y0 + j * dy;

            int index = j * width + i;
            output[index] = mandel(x, y, maxIterations);
    }
}

但是通过prog2我们知道,SIMD的实现方式总是会带着一个向量mask(这里是8位),mask各个bit的值越是不趋同,越是产生更多的execution divergence。如下图所示,那些在黑白相间的区间执行的SIMD计算,会产生更多的效率损失,这就是加速比达不到理论值的原因。可以对比view2,它黑白相间的区间比view1更多,所以它的加速比更低。

image-20260625232657265

part2

使用ispctask的形式完成计算,启用多个task相当于让ispc利用了多个核。我的实验环境是4核8线程,按照part1 view1的加速比4.86作为参考,理论上当task数量等于8时,我们就可达到最大加速比 4.86 * 8 = 38.88

试验结果如下表所示,随着tasknum的增大,加速比逐渐增大,当tasknum达到16时,达到最大加速比只有35.78,随后加速比随着tasknum的增加而减少。

tasknum task_ispc (ms) ispc (ms) serial (ms) ispc speedup task speedup
2 24.061 47.954 231.372 4.824874 9.616059
4 19.038 47.006 233.221 4.961516 12.250289
8 11.902 48.002 231.746 4.827841 19.471181
16 7.736 48.183 230.771 4.789469 29.830791
32 6.499 47.856 232.542 4.859203 35.781197
80 7.314 47.435 230.633 4.862085 31.533087
160 7.477 47.859 230.611 4.818550 30.842718

解释如下:加速比最大值没有达到预期是因为负载不平衡(prog1已经得到该结论);而tasknum为8时没有达到最大加速比,反而在16时才达到,为什么?因为更多的线程数使得负载分配更加均衡,但是当线程数量过多时,其带来的开销会大于收益,所以加速比会渐渐减少。

Prog4 Iterative sqrt

使用ispc实现牛顿迭代法解平方根,迭代次数与值的关系入下:

image-20260626000058048

主要回答Readme中的三个问题:

Q1:Build and run sqrt. Report the ISPC implementation speedup for single CPU core (no tasks) and when using all cores (with tasks). What is the speedup due to SIMD parallelization? What is the speedup due to multi-core parallelization?

实验结果如下:

~~~md
[sqrt serial]:          [778.245] ms
[sqrt ispc]:            [192.859] ms
[sqrt task ispc]:       [25.875] ms
                                (4.04x speedup from ISPC)
                                (30.08x speedup from task ISPC)
~~~



单核CPU到达4.04加速比,多核CPU达到30.08加速比。4.04是SIMD并行的共性,(30.08 / 4.04 =) 7.44是多核并行的贡献

Q2:Modify the contents of the array values to improve the relative speedup of the ISPC implementations. Construct a specifc input that maximizes speedup over the sequential version of the code and report the resulting speedup achieved (for both the with- and without-tasks ISPC implementations). Does your modification improve SIMD speedup? Does it improve multi-core speedup (i.e., the benefit of moving from ISPC without-tasks to ISPC with tasks)? Please explain why.

为达到最大加速比,应当尽可能让所有运算通道实现负载均衡,同时负载量要足够大,以此抵消内存读取、上下文切换等和计算本身无关的开销

因此,最优取值应无限接近于 3。但受限于浮点型数据的精度,我最终选用了2.999999f

for (unsigned int i=0; i<N; i++)
{
    values[i] = 2.999999f;
}

测试结果如下:

[sqrt serial]:          [4595.527] ms
[sqrt ispc]:            [701.672] ms
[sqrt task ispc]:       [90.519] ms
                                (6.55x speedup from ISPC)
                                (50.77x speedup from task ISPC)

也就是说,SIMD 单指令多数据带来 6.55 倍加速,多核并行额外带来 7.75 倍加速。两种并行的加速比相对Q1来说都提升了,SIMD提升是因为负载均衡,多核提升是因为计算量大。

Q3 :Construct a specific input for sqrt that minimizes speedup for ISPC (without-tasks) over the sequential version of the code. Describe this input, describe why you chose it, and report the resulting relative performance of the ISPC implementations. What is the reason for the loss in efficiency? (keep in mind we are using the --target=avx2 option for ISPC, which generates 8-wide SIMD instructions).

若要把加速比降到最低,就需要在单次向量运算中引入execution divergence。同时我知道:当计算负载(相对于内存访问以及其他操作)越小时,最终的加速效果就越差。

因此在每 8 个元素的计算任务里,将其中 1 个元素赋值为2.999999f,剩余 7 个元素赋值为1.f,这样配置能够得到最低的加速比:

for (unsigned int i=0; i<N; i++)
{
    if (i % 8 == 0)
    {
        values[i] = 2.999999f;
    }
    else 
    {
        values[i] = 1.f;
    }
}

最终结果如下:

[sqrt serial]:          [566.620] ms
[sqrt ispc]:            [675.777] ms
[sqrt task ispc]:       [90.690] ms
                                (0.84x speedup from ISPC)
                                (6.25x speedup from task ISPC)

Prog5 BLAS saxpy

结果如下,开启 ISPC 多任务并行后仅能获得 1.30 倍的加速比,我认为主要原因是该任务属于 IO 密集型负载

[saxpy ispc]:           [11.750] ms     [25.363] GB/s   [3.404] GFLOPS
[saxpy task ispc]:      [9.065] ms      [32.875] GB/s   [4.412] GFLOPS
                                (1.30x speedup from use of tasks)

理论计算分析:我的 CPU 型号为:11 代英特尔酷睿 i5-11400H,主频 2.70GHz,且 WSL 将可用核心数限制为 4 核。

理论峰值浮点算力 = 核心数 × 主频 (GHz) × 单周期浮点运算量,即:4 × 2.7 × 8(实际值可能更高)= 86.4 GFLOPS。

要跑满这一峰值算力,所需内存带宽应为:86.4 GFPLOPS * 1FPLOPS * 4BYTES(32FP) = 345.6GB/s

而该 官方给出该CPU的理论内存带宽仅为 51.2GB / s,由此可见,现代 CPU 的运算能力并不是性能瓶颈,真正的瓶颈在于内存带宽。

Prog6 Making K-Means Faster

该实验需要优化kmeans算法,readme先让我们测程序的各个部分耗时,然后选择一个进行优化

初始测试结果

computeAssignments cost 341.352515 ms
computeCentroids cost 61.992996 ms
computeCost cost 103.041950 ms
[Total Time]: 10900.018 ms

我选择优化 computeAssignments

void computeAssignments(WorkerArgs *const args) {
  double *minDist = new double[args->M];
  // Initialize arrays
  for (int m =0; m < args->M; m++) {
    minDist[m] = 1e30;
    args->clusterAssignments[m] = -1;
  }
  int threadNum = 7;
  std::thread workers[threadNum];
  int step = (args->M + 1) / (threadNum + 1);

  for (int i = 0; i < threadNum; i++)
  {
    workers[i] = std::thread(threadWork, args, step * i, step * (i + 1), minDist);
  }
  // Assign datapoints to closest centroids
  threadWork(args, step * threadNum, args->M, minDist);
  for (int i = 0; i < threadNum; i++)
  {
    workers[i].join();
  }
  delete[] minDist;
}

优化后:

computeAssignments cost 39.449921 ms 
computeCentroids cost 55.315588 ms 
computeCost cost 92.758145 ms 
[Total Time]: 4291.821 ms

最没意思的一集

Assignment 2

实现一个任务执行的库(使用多核多线程),api类似ISPC。分为两个部分,第一部分的任务无依赖关系,且API是同步的,第二部分则需要维护任务的依赖关系,API是异步的

Part A Synchronous Bulk Task Launch

实验知道引导我们依次实现三种任务处理类:

  • TaskSystemParallelSpawn,无线程池
  • TaskSystemParallelThreadPoolSpinning,有线程池,线程在没有任务时自旋
  • TaskSystemParallelThreadPoolSleeping,有线程池,线程在没有任务时睡眠
TaskSystemParallelSpawn

无线程池的实现,调用 TaskSystemParallelSpawn::run 时创建线程,执行完任务后再销毁。

总体行为较简单,注意这里也要尽可能在任务颗粒度及其带来的开销间进行平衡

void TaskSystemParallelSpawn::run(IRunnable* runnable, int num_total_tasks) {
    int multipler = 16;
    int num_slice = m_num_threads * multipler;
    int num_tasks_each_slice = num_total_tasks / (num_slice);
    if (num_tasks_each_slice <= 0)
    {
        num_tasks_each_slice = 1;
    }
    for (int i = 0; i < m_num_threads - 1; i++) {
        m_threads[i] = std::thread(TaskSystemParallelSpawnThreadWorker,i, runnable, num_total_tasks, i * num_tasks_each_slice, num_tasks_each_slice, m_num_threads * num_tasks_each_slice);
    }
    TaskSystemParallelSpawnThreadWorker(m_num_threads - 1, runnable, num_total_tasks, (m_num_threads - 1) * num_tasks_each_slice, num_tasks_each_slice, m_num_threads * num_tasks_each_slice);
    if (num_slice * num_tasks_each_slice < num_total_tasks)
    {
        for (int i = num_slice * num_tasks_each_slice; i < num_total_tasks; i++)
        {
            runnable->runTask(i, num_total_tasks);
        }
    }

    for (int i = 0; i < m_num_threads - 1; i++) {
        m_threads[i].join();
    }
}

测试结果如下:

                                        STUDENT   REFERENCE   PERF?
[Serial]                                1630.528  1609.641    1.01  (OK)
[Parallel + Always Spawn]               403.927   360.977     1.12  (OK)
[Parallel + Thread Pool + Spin]         1519.358  346.356     4.39  (NOT OK)
[Parallel + Thread Pool + Sleep]        1607.045  342.001     4.70  (NOT OK)
// ......
Overall performance results
[Serial]                                : All passed Perf
[Parallel + Always Spawn]               : All passed Perf
[Parallel + Thread Pool + Spin]         : Perf did not pass all tests
[Parallel + Thread Pool + Sleep]        : Perf did not pass all tests
TaskSystemParallelThreadPoolSpinning

在TaskSystemParallelThreadPoolSpinning初始化时创建线程池,在其销毁时再join所有的线程,当线程没有任务可做时,自旋

我为每个线程配了一个WorkQueue,WorkQueue拥有如下API:

class WorkQueue
{
public:
    WorkQueue() = default;
    void InQueueTask(TaskDescription& taskdes);
    bool DeQueueTask(TaskDescription& out);
    TaskDescription& Front();
    bool IsEmpty();
private:
    std::queue<TaskDescription> m_queue;
    std::mutex m_lock;
};
class TaskDescription
{
public:
    int assign_from;
    int assign_to;
    IRunnable* runnable;
    int num_total_tasks;
    TaskDescription() : assign_from(0), assign_to(0), runnable(nullptr), num_total_tasks(0) {}
    TaskDescription(int from, int to, IRunnable* r, int total): 
            assign_from(from), assign_to(to), runnable(r), num_total_tasks(total) {}
};

TaskSystemParallelThreadPoolSpinning在初始化时就创建一个线程池,销毁时再join:

TaskSystemParallelThreadPoolSpinning::TaskSystemParallelThreadPoolSpinning(int num_threads): 
    ITaskSystem(num_threads), 
    m_taskState(),
    m_workQuque(num_threads),
    m_threads(num_threads - 1), 
    m_num_threads(num_threads)
{
    m_taskState.m_killed.store(false, std::memory_order_release);
    m_taskState.m_RemainingTask.store(0, std::memory_order_release);
    for (size_t i = 0; i < m_threads.size(); i++)
    {
        m_threads[i] = std::thread(TaskSystemParallelThreadPoolSpinningWorker, static_cast<int>(i), this);
    }
}
TaskSystemParallelThreadPoolSpinning::~TaskSystemParallelThreadPoolSpinning() {
    // 给自旋线程一个kill信号
    m_taskState.m_killed.store(true, std::memory_order_release);
    for (int i = 0; i < m_num_threads - 1; i++) {
        m_threads[i].join();
    }
}

每个线程循环地调用DoWork,该函数检查自己的WorkQueue是否有任务,如果没有从其他线程的WorkQueue中“偷”一个。注意这种任务分配方式属于动态分配,如果还使用静态分配,最后的测试结果有好多个不达标,可以参考代码仓库中的ANSWER.md,这里记录了多种实现方式及其测试结果。

void TaskSystemParallelThreadPoolSpinningWorker(int index, TaskSystemParallelThreadPoolSpinning* taskSystem)
{
    // 每个线程循环调用DoWork()
    while (true)
    {
        // 如果有kill信号,结束循环退出
        if (taskSystem->m_taskState.m_killed.load(std::memory_order_acquire))
        {
            break;
        } 
        DoWork(index, taskSystem);
    }
}

void DoWork(int index, TaskSystemParallelThreadPoolSpinning* taskSystem)
{
    TaskDescription taskDes;
    // 先从自己的workqueue中找任务,如果找不到再从别的WorkQueue中找
    bool doSelf = taskSystem->m_workQuque[index].DeQueueTask(taskDes);
    if (!doSelf)
    {
        // steal from others
        for (int stolenIdx = (index + 1) % (taskSystem->m_workQuque.size()); 
            stolenIdx != index; 
            stolenIdx = (stolenIdx + 1) % (taskSystem->m_workQuque.size()))
        {
            if (taskSystem->m_workQuque[stolenIdx].DeQueueTask(taskDes))
            {
                break;
            }
        }
    }

    if (taskDes.runnable)
    {
        for (int i = taskDes.assign_from; i < taskDes.assign_to; i++)
        {
            taskDes.runnable->runTask(i, taskDes.num_total_tasks);
        }
        int doneNum = taskDes.assign_to - taskDes.assign_from;
        taskSystem->m_taskState.m_RemainingTask.fetch_sub(doneNum, std::memory_order_seq_cst);
    }
}

最后,TaskSystemParallelThreadPoolSpinningWorker的调用接口如下,注意由于API是同步的,所以调用run的线程也可以不停地获取任务并执行任务,直到任务全部执行完毕

void TaskSystemParallelThreadPoolSpinning::run(IRunnable* runnable, int num_total_tasks) 
{
	// 1. 切分任务,为每个Workqueue填充任务描述符
    int multipler = 8;
    int num_slice = m_num_threads * multipler;
    int num_tasks_each_slice = num_total_tasks / (num_slice);
    if (num_tasks_each_slice <= 0)
    {
        num_tasks_each_slice = 1;
    }
    m_taskState.m_RemainingTask.store(num_total_tasks, std::memory_order_release);
    int assignedIdx = 0; 
    for (int i = 0; i < num_total_tasks; i += num_tasks_each_slice)
    {
        TaskDescription taskDes{i, i + num_tasks_each_slice, runnable, num_total_tasks};
        if (i + num_tasks_each_slice > num_total_tasks)
        {
            taskDes.assign_to = num_total_tasks;
        }
        m_workQuque[assignedIdx].InQueueTask(taskDes);
        assignedIdx = (assignedIdx + 1) % m_workQuque.size();
    }
	// 2. 分配好任务后,主线程也循环调用DoWork执行运算知道系统的所有任务执行完毕
    while (m_taskState.m_RemainingTask.load(std::memory_order_acquire) > 0)
    {
        DoWork(m_num_threads - 1, this);
    }
    
}

测试结果如下,详细结果见ANSWER.md, 不知道为什么Parallel + Always Spawn会有回退

Overall performance results
[Serial]                                : All passed Perf
[Parallel + Always Spawn]               : Perf did not pass all tests
[Parallel + Thread Pool + Spin]         : All passed Perf
[Parallel + Thread Pool + Sleep]        : Perf did not pass all tests
TaskSystemParallelThreadPoolSleeping

工作线程在没有任务时不再自旋,而是进入睡眠,这需要引入多个同步变量,为了方便阅读,将他们聚合到一个类TaskState里

class TaskSystemParallelThreadPoolSleeping: public ITaskSystem {
    public:
        class TaskState
        {
            public:
            std::atomic<bool> m_killed;  // 整个系统是否退出的标志位
            std::atomic<int> m_RemainingTask; // 整个系统还有多少任务
            // 是否还有任务的条件变量 + 锁,用于线程的睡眠和唤醒
            std::mutex  m_hasTaksLk;
            std::condition_variable m_hasTaksCv; 
        };
    // ...
            TaskState m_taskState;
        std::vector<WorkQueue> m_workQuque;

    private:
        std::vector<std::thread> m_threads; // m_threads.size() + 1 = m_num_threads;
        int m_num_threads;
}

同上一个实现,系统初始化和销毁时分别创建和销毁线程池:

TaskSystemParallelThreadPoolSleeping::TaskSystemParallelThreadPoolSleeping(int num_threads):
    ITaskSystem(num_threads), 
    m_taskState(),
    m_workQuque(num_threads),
    m_threads(num_threads - 1), 
    m_num_threads(num_threads) 
{
    m_taskState.m_killed.store(false, std::memory_order_release);
    m_taskState.m_RemainingTask.store(0, std::memory_order_release);
    for (size_t i = 0; i < m_threads.size(); i++)
    {
        m_threads[i] = std::thread(TaskSystemParallelThreadPoolSleepingWorker, static_cast<int>(i), this);
    }
}

TaskSystemParallelThreadPoolSleeping::~TaskSystemParallelThreadPoolSleeping() {
    m_taskState.m_killed.store(true, std::memory_order_release);
    std::unique_lock<std::mutex> lc(m_taskState.m_hasTaksLk);
    m_taskState.m_hasTaksCv.notify_all(); // 这里的唤醒是必要的
    lc.unlock();
    for (int i = 0; i < m_num_threads - 1; i++) {
        m_threads[i].join();
    }
}


其他函数的实现就不一一例举了, 最需要注意点就是线程之间的 睡眠唤醒 竞争,其他的和自旋系统差不多。试验结果如下:

Overall performance results
[Serial]                                : All passed Perf
[Parallel + Always Spawn]               : All passed Perf
[Parallel + Thread Pool + Spin]         : All passed Perf
[Parallel + Thread Pool + Sleep]        : All passed Perf

PartB Supporting Execution of Task Graphs

相比于PartA的要求,PartB的任务有依赖关系,且API是异步的:

  1. TaskSystemParallelThreadPoolSleeping::runAsyncWithDeps :往系统里赛任务,可连续调用
  2. TaskSystemParallelThreadPoolSleeping::sync:塞完所有任务后被调用,等待系统的所有任务完成后返回

复杂度一下子上来了,首先系统API是异步的,那么我需要两个queue,一个是workqueue一个是waitingqueue,外部系统调用时runAsyncWithDeps ,把它塞到在waitingqueue中,然后立刻返回。工作线程发现workqueue中没有任务时,再从waitingqueue中取出任务放到workqueue中。

此外,维护依赖关系,还需要一个数据结构去记录依赖,以便于在任务结束时更新其他任务的依赖数量。

如下所示,是本系统用到的所有数据结构,

class TaskSystemParallelThreadPoolSleeping: public ITaskSystem {
    public:
        class TaskState
        {
            public:
                std::atomic<int> m_RemainingTask;
                std::atomic<bool> m_killed;

                std::mutex  m_hasTaksLk;
                std::condition_variable m_hasTaksCv;
				// 主要是给sync用的
                std::mutex m_allTasksDoneLk;
                std::condition_variable m_allTasksDoneCv;
        };
        TaskSystemParallelThreadPoolSleeping(int num_threads);
        ~TaskSystemParallelThreadPoolSleeping();
        const char* name();
        void run(IRunnable* runnable, int num_total_tasks);
        TaskID runAsyncWithDeps(IRunnable* runnable, int num_total_tasks,
                                const std::vector<TaskID>& deps);
        void SliceAndMoveTheTaskToWorkQueue(const TaskDesc& taskDes);
        bool GetWorkFromWaitingQueue(int index);
        void sync();

        void TaskSystemParallelThreadPoolSleepingWorker(int index);


        TaskID m_nextTaskId; // 相当于uniqueid,为每个任务标号
        std::vector<std::thread> m_threads;
        int m_num_threads;

        TaskState m_taskState;
   
        std::vector<WorkQueue>  m_workQuque; // task slices that are ready to exec and are distributed to m_num_threads
    
        std::vector<TaskDesc*>   m_waitingQueue;
        std::mutex              m_waitingQueueLck;

        std::mutex              m_taskRecordLck; //  记录依赖关系
        std::unordered_map<TaskID, std::unique_ptr<TaskDesc>> m_taskRecords;

};

WorkQueue的API与上一版本是实现类似:

class WorkQueue
{
public:
    WorkQueue() = default;
    void InQueueTask(TaskSliceDesc& taskdes);
    bool DeQueueTask(TaskSliceDesc& out);
    TaskSliceDesc& Front();
    bool IsEmpty();
private:
    std::queue<TaskSliceDesc> m_queue;
    std::mutex m_lock;
};

TaskDesc直接描述runAsyncWithDeps传递过来的外部任务,包括它的总任务数,依赖关系等,此外为了方便维护,系统内部还添加了 task_not_done_num、

deped_has_not_been_done_num、sufs(记录依赖这个任务的任务)等字段

class TaskDesc
{
    public:
        IRunnable* runnable;
        int num_total_tasks;
        TaskID taskId;

        std::atomic<int> task_not_done_num; // decreased from num_total_tasks
        std::atomic<int> deped_has_not_been_done_num; // Only when all dependent tasks are completed can this task moved into workqueue
        std::vector<TaskID> deps;
        std::set<TaskID> sufs; // Tasks that depend on this task
        std::mutex taskDesLck;
    TaskDesc(): runnable(nullptr), num_total_tasks(0), taskId(-1), task_not_done_num(-1) {}
    TaskDesc(IRunnable* r, int num_total_tasks, TaskID taskId, const std::vector<TaskID>& indeps):
            runnable(r), num_total_tasks(num_total_tasks), taskId(taskId), 
            task_not_done_num(num_total_tasks), deped_has_not_been_done_num(0), deps(indeps),  sufs(), taskDesLck() {}
};


TaskSliceDesc是TaskDesc被进行切分后的存在形式,WorkQueue直接存储TaskSliceDesc而不是TaskDesc,好处就是工作线程拿过来就用,不需要再进行切分操作。

class TaskSliceDesc
{
    public:
        int assign_from;
        int assign_to;
        IRunnable* runnable;
        int num_total_tasks;
        TaskID taskid_belongs_to;
        TaskSliceDesc() : assign_from(0), assign_to(0), runnable(nullptr), num_total_tasks(0), taskid_belongs_to(-1) {}
        TaskSliceDesc(int from, int to, IRunnable* r, int total, TaskID taskIdBelongsTo): 
                assign_from(from), assign_to(to), runnable(r), num_total_tasks(total), taskid_belongs_to(taskIdBelongsTo) {}
        
};

具体的函数实现,这里就讲解下工作线程的实现吧,自认为整个partB的实现是比较复杂的(和其他人的参考作业相比相当复杂。因为我看了Lecture5后,就先入为主地认为这里需要用workqueue实现,但其实不用也能实现)就不详细介绍了。

线程工作函数如下所示,它依然是一个while循环,他主要做四件事情:

  1. 判断当前系统是否销毁,若是则退出循环
  2. 调用DoWork函数进行真正的任务运行,如果其返回值非零,表示WorkQueue中还有任务并已经执行完成,则continue进入下个循环,若返回值是0则表示所有线程的WorkQueue中都没有任务了,进入第三步骤
  3. 调用GetWorkFromWaitingQueue函数,将任务从waitingqueue搬到workqueue,并且做完所有的任务切分和任务分配工作,做完后进入下一个循环。若其返回值为0,表示waitingqueue和workqueue都为空,进入步骤4
  4. 等待其他任务下发到系统,注意这里睡眠唤醒竞争,这是再判断一次m_killed信号才调用conditionwait的原因
void TaskSystemParallelThreadPoolSleeping::TaskSystemParallelThreadPoolSleepingWorker(int index)
{
    while (true)
    {
        if (m_taskState.m_killed.load(std::memory_order_acquire))
        {
            break;
        } 

        if (DoWork(index, this))
        {
            continue;
        }
        else if (GetWorkFromWaitingQueue(index))
        {
            continue;
        }
        else
        {
            std::unique_lock<std::mutex> lc(m_taskState.m_hasTaksLk);
            // avoid sleep/notify condition when destroying the thread system
            if (!m_taskState.m_killed.load(std::memory_order_acquire) && m_taskState.m_RemainingTask.load() == 0)
            {
                debugprint("threadidx %d wait for work\n", index);
                m_taskState.m_hasTaksCv.wait(lc);
            }
        }
    }
}

测试结果如下所示,super_light这个实验没有达标。原因是这个任务非常轻量,但我的实现设计大线程间通信、同步、锁的操作,所以在计算量小的任务中没有取得很好的结果

runtasks_ref
Linux x86_64
================================================================================
Running task system grading harness... (22 total tests)
  - Detected CPU with 8 execution contexts
  - Task system configured to use at most 8 threads
================================================================================
================================================================================
Executing test: super_super_light...
Reference binary: ./runtasks_ref_linux
Results for: super_super_light
                                        STUDENT   REFERENCE   PERF?
[Serial]                                4.59      5.176       0.89  (OK)
[Parallel + Always Spawn]               5.513     120.967     0.05  (OK)
[Parallel + Thread Pool + Spin]         5.387     23.058      0.23  (OK)
[Parallel + Thread Pool + Sleep]        47.596    39.051      1.22  (OK)
================================================================================
Executing test: super_super_light_async...
Reference binary: ./runtasks_ref_linux
Results for: super_super_light_async
                                        STUDENT   REFERENCE   PERF?
[Serial]                                4.623     5.193       0.89  (OK)
[Parallel + Always Spawn]               5.351     122.087     0.04  (OK)
[Parallel + Thread Pool + Spin]         5.314     18.023      0.29  (OK)
[Parallel + Thread Pool + Sleep]        11.853    24.684      0.48  (OK)
================================================================================
Executing test: super_light...
Reference binary: ./runtasks_ref_linux
Results for: super_light
                                        STUDENT   REFERENCE   PERF?
[Serial]                                63.879    66.577      0.96  (OK)
[Parallel + Always Spawn]               63.511    129.195     0.49  (OK)
[Parallel + Thread Pool + Spin]         62.504    30.364      2.06  (NOT OK)
[Parallel + Thread Pool + Sleep]        67.457    41.781      1.61  (NOT OK)
================================================================================
Executing test: super_light_async...
Reference binary: ./runtasks_ref_linux
Results for: super_light_async
                                        STUDENT   REFERENCE   PERF?
[Serial]                                63.088    67.028      0.94  (OK)
[Parallel + Always Spawn]               62.572    127.858     0.49  (OK)
[Parallel + Thread Pool + Spin]         62.213    25.487      2.44  (NOT OK)
[Parallel + Thread Pool + Sleep]        29.413    34.094      0.86  (OK)
================================================================================
Executing test: ping_pong_equal...
Reference binary: ./runtasks_ref_linux
Results for: ping_pong_equal
                                        STUDENT   REFERENCE   PERF?
[Serial]                                1045.914  1104.507    0.95  (OK)
[Parallel + Always Spawn]               1044.569  312.782     3.34  (NOT OK)
[Parallel + Thread Pool + Spin]         1040.416  245.787     4.23  (NOT OK)
[Parallel + Thread Pool + Sleep]        253.139   255.036     0.99  (OK)
================================================================================
Executing test: ping_pong_equal_async...
Reference binary: ./runtasks_ref_linux
Results for: ping_pong_equal_async
                                        STUDENT   REFERENCE   PERF?
[Serial]                                1034.08   1110.039    0.93  (OK)
[Parallel + Always Spawn]               1033.909  316.404     3.27  (NOT OK)
[Parallel + Thread Pool + Spin]         1025.556  255.766     4.01  (NOT OK)
[Parallel + Thread Pool + Sleep]        217.495   242.501     0.90  (OK)
================================================================================
Executing test: ping_pong_unequal...
Reference binary: ./runtasks_ref_linux
Results for: ping_pong_unequal
                                        STUDENT   REFERENCE   PERF?
[Serial]                                1555.376  1565.179    0.99  (OK)
[Parallel + Always Spawn]               1559.839  393.022     3.97  (NOT OK)
[Parallel + Thread Pool + Spin]         1550.855  331.264     4.68  (NOT OK)
[Parallel + Thread Pool + Sleep]        334.977   333.978     1.00  (OK)
================================================================================
Executing test: ping_pong_unequal_async...
Reference binary: ./runtasks_ref_linux
Results for: ping_pong_unequal_async
                                        STUDENT   REFERENCE   PERF?
[Serial]                                1556.645  1567.986    0.99  (OK)
[Parallel + Always Spawn]               1556.532  397.567     3.92  (NOT OK)
[Parallel + Thread Pool + Spin]         1551.312  331.157     4.68  (NOT OK)
[Parallel + Thread Pool + Sleep]        297.209   301.373     0.99  (OK)
================================================================================
Executing test: recursive_fibonacci...
Reference binary: ./runtasks_ref_linux
Results for: recursive_fibonacci
                                        STUDENT   REFERENCE   PERF?
[Serial]                                782.421   1311.051    0.60  (OK)
[Parallel + Always Spawn]               780.507   240.921     3.24  (NOT OK)
[Parallel + Thread Pool + Spin]         780.502   250.949     3.11  (NOT OK)
[Parallel + Thread Pool + Sleep]        158.599   227.56      0.70  (OK)
================================================================================
Executing test: recursive_fibonacci_async...
Reference binary: ./runtasks_ref_linux
Results for: recursive_fibonacci_async
                                        STUDENT   REFERENCE   PERF?
[Serial]                                785.134   1308.337    0.60  (OK)
[Parallel + Always Spawn]               782.655   238.588     3.28  (NOT OK)
[Parallel + Thread Pool + Spin]         779.706   234.047     3.33  (NOT OK)
[Parallel + Thread Pool + Sleep]        149.887   220.149     0.68  (OK)
================================================================================
Executing test: math_operations_in_tight_for_loop...
Reference binary: ./runtasks_ref_linux
Results for: math_operations_in_tight_for_loop
                                        STUDENT   REFERENCE   PERF?
[Serial]                                492.855   518.46      0.95  (OK)
[Parallel + Always Spawn]               499.716   514.302     0.97  (OK)
[Parallel + Thread Pool + Spin]         502.062   166.016     3.02  (NOT OK)
[Parallel + Thread Pool + Sleep]        387.519   260.237     1.49  (OK)
================================================================================
Executing test: math_operations_in_tight_for_loop_async...
Reference binary: ./runtasks_ref_linux
Results for: math_operations_in_tight_for_loop_async
                                        STUDENT   REFERENCE   PERF?
[Serial]                                500.134   506.514     0.99  (OK)
[Parallel + Always Spawn]               502.658   498.004     1.01  (OK)
[Parallel + Thread Pool + Spin]         497.772   155.128     3.21  (NOT OK)
[Parallel + Thread Pool + Sleep]        250.697   216.439     1.16  (OK)
================================================================================
Executing test: math_operations_in_tight_for_loop_fewer_tasks...
Reference binary: ./runtasks_ref_linux
Results for: math_operations_in_tight_for_loop_fewer_tasks
                                        STUDENT   REFERENCE   PERF?
[Serial]                                508.337   514.144     0.99  (OK)
[Parallel + Always Spawn]               511.246   484.791     1.05  (OK)
[Parallel + Thread Pool + Spin]         505.98    174.051     2.91  (NOT OK)
[Parallel + Thread Pool + Sleep]        416.984   288.946     1.44  (OK)
================================================================================
Executing test: math_operations_in_tight_for_loop_fewer_tasks_async...
Reference binary: ./runtasks_ref_linux
Results for: math_operations_in_tight_for_loop_fewer_tasks_async
                                        STUDENT   REFERENCE   PERF?
[Serial]                                509.615   523.233     0.97  (OK)
[Parallel + Always Spawn]               503.328   488.903     1.03  (OK)
[Parallel + Thread Pool + Spin]         507.326   99.029      5.12  (NOT OK)
[Parallel + Thread Pool + Sleep]        90.798    91.546      0.99  (OK)
================================================================================
Executing test: math_operations_in_tight_for_loop_fan_in...
Reference binary: ./runtasks_ref_linux
Results for: math_operations_in_tight_for_loop_fan_in
                                        STUDENT   REFERENCE   PERF?
[Serial]                                264.016   268.041     0.98  (OK)
[Parallel + Always Spawn]               262.338   94.796      2.77  (NOT OK)
[Parallel + Thread Pool + Spin]         263.522   66.412      3.97  (NOT OK)
[Parallel + Thread Pool + Sleep]        85.942    73.893      1.16  (OK)
================================================================================
Executing test: math_operations_in_tight_for_loop_fan_in_async...
Reference binary: ./runtasks_ref_linux
Results for: math_operations_in_tight_for_loop_fan_in_async
                                        STUDENT   REFERENCE   PERF?
[Serial]                                263.967   269.124     0.98  (OK)
[Parallel + Always Spawn]               261.748   93.892      2.79  (NOT OK)
[Parallel + Thread Pool + Spin]         260.641   52.182      4.99  (NOT OK)
[Parallel + Thread Pool + Sleep]        47.044    49.104      0.96  (OK)
================================================================================
Executing test: math_operations_in_tight_for_loop_reduction_tree...
Reference binary: ./runtasks_ref_linux
Results for: math_operations_in_tight_for_loop_reduction_tree
                                        STUDENT   REFERENCE   PERF?
[Serial]                                262.318   267.759     0.98  (OK)
[Parallel + Always Spawn]               258.712   60.455      4.28  (NOT OK)
[Parallel + Thread Pool + Spin]         258.588   56.979      4.54  (NOT OK)
[Parallel + Thread Pool + Sleep]        54.306    56.539      0.96  (OK)
================================================================================
Executing test: math_operations_in_tight_for_loop_reduction_tree_async...
Reference binary: ./runtasks_ref_linux
Results for: math_operations_in_tight_for_loop_reduction_tree_async
                                        STUDENT   REFERENCE   PERF?
[Serial]                                258.274   262.23      0.98  (OK)
[Parallel + Always Spawn]               257.322   59.635      4.31  (NOT OK)
[Parallel + Thread Pool + Spin]         256.878   49.653      5.17  (NOT OK)
[Parallel + Thread Pool + Sleep]        44.615    45.9        0.97  (OK)
================================================================================
Executing test: spin_between_run_calls...
Reference binary: ./runtasks_ref_linux
Results for: spin_between_run_calls
                                        STUDENT   REFERENCE   PERF?
[Serial]                                277.426   480.404     0.58  (OK)
[Parallel + Always Spawn]               279.996   244.731     1.14  (OK)
[Parallel + Thread Pool + Spin]         280.246   283.766     0.99  (OK)
[Parallel + Thread Pool + Sleep]        178.431   246.618     0.72  (OK)
================================================================================
Executing test: spin_between_run_calls_async...
Reference binary: ./runtasks_ref_linux
Results for: spin_between_run_calls_async
                                        STUDENT   REFERENCE   PERF?
[Serial]                                279.213   476.81      0.59  (OK)
[Parallel + Always Spawn]               279.879   244.729     1.14  (OK)
[Parallel + Thread Pool + Spin]         281.76    279.576     1.01  (OK)
[Parallel + Thread Pool + Sleep]        180.042   245.616     0.73  (OK)
================================================================================
Executing test: mandelbrot_chunked...
Reference binary: ./runtasks_ref_linux
Results for: mandelbrot_chunked
                                        STUDENT   REFERENCE   PERF?
[Serial]                                418.071   414.412     1.01  (OK)
[Parallel + Always Spawn]               420.103   57.948      7.25  (NOT OK)
[Parallel + Thread Pool + Spin]         415.18    63.548      6.53  (NOT OK)
[Parallel + Thread Pool + Sleep]        57.813    58.356      0.99  (OK)
================================================================================
Executing test: mandelbrot_chunked_async...
Reference binary: ./runtasks_ref_linux
Results for: mandelbrot_chunked_async
                                        STUDENT   REFERENCE   PERF?
[Serial]                                414.91    417.726     0.99  (OK)
[Parallel + Always Spawn]               416.311   57.788      7.20  (NOT OK)
[Parallel + Thread Pool + Spin]         419.439   61.685      6.80  (NOT OK)
[Parallel + Thread Pool + Sleep]        57.561    58.114      0.99  (OK)
================================================================================
Overall performance results
[Serial]                                : All passed Perf
[Parallel + Always Spawn]               : Perf did not pass all tests
[Parallel + Thread Pool + Spin]         : Perf did not pass all tests
[Parallel + Thread Pool + Sleep]        : Perf did not pass all tests

Assignment 3

使用Cuda进行加速运算的一个project,由于不熟悉Cuda,所以写起来比较吃力。

环境:

此外,ref程序需要特定的cudaruntime版本,但是太高的版本又和我本地环境不适配,所以checkout了23年的某个commit,然后一把copy了过来。

SAXPY

练手的,用Cuda再实现一遍asst1中的SAXPY,熟悉cudaMalloc, cudaMemcpy, cudaFree, cudaDeviceSynchronize 以及如何launch一个cuda核函数。

试验结果对比asst1中的ispc如下所示,需要将asst1中的数据量改成本实验的数据量对比:

ispc实现:

[saxpy serial]:         [76.086] ms     [19.585] GB/s   [2.629] GFLOPS
[saxpy ispc]:           [75.910] ms     [19.630] GB/s   [2.635] GFLOPS
[saxpy task ispc]:      [57.970] ms     [25.705] GB/s   [3.450] GFLOPS

cuda实现:

Found 1 CUDA devices
Device 0: NVIDIA GeForce RTX 3050 Laptop GPU
   SMs:        16
   Global mem: 4096 MB
   CUDA Cap:   8.6
---------------------------------------------------------
Running 3 timing tests:
Effective BW by CUDA saxpy: 233.917 ms          [4.778 GB/s]
Effective BW of kernel calc: 10.398 ms          [107.480 GB/s]
Effective BW by CUDA saxpy: 173.751 ms          [6.432 GB/s]
Effective BW of kernel calc: 6.959 ms           [160.595 GB/s]
Effective BW by CUDA saxpy: 192.386 ms          [5.809 GB/s]
Effective BW of kernel calc: 7.489 ms           [149.233 GB/s]

对比计算时间,cuda的计算效率高很多

[saxpy task ispc]:      [57.970] ms     [25.705] GB/s   [3.450] GFLOPS
Effective BW of kernel calc: 7.489 ms           [149.233 GB/s]

但是cuda内存传输的时间非常耗时,贷款达149.233 GB/s 比较接近理论带宽(193GB/s),说明程序瓶颈在带宽上。

Scan

实现exclusive_scan,具体实现可以参考视频教程Lecture8,PA的readme也有些对应的伪代码:

void exclusive_scan_iterative(int* start, int* end, int* output) {
    int N = end - start;
    memmove(output, start, N*sizeof(int));
    // upsweep phase
    for (int two_d = 1; two_d <= N/2; two_d*=2) {
        int two_dplus1 = 2*two_d;
        parallel_for (int i = 0; i < N; i += two_dplus1) {
            output[i+two_dplus1-1] += output[i+two_d-1];
        }
    }
    output[N-1] = 0;
    // downsweep phase
    for (int two_d = N/2; two_d >= 1; two_d /= 2) {
        int two_dplus1 = 2*two_d;
        parallel_for (int i = 0; i < N; i += two_dplus1) {
            int t = output[i+two_d-1];
            output[i+two_d-1] = output[i+two_dplus1-1];
            output[i+two_dplus1-1] += t;
        }
    }
}

不是很难,照着实现就行,只有一个注意点就是注意爆int,我踩过坑所以把相应的kernel函数中int全改成了long long

__global__ void
upsweep_kernel(int N, int* inputArray, int stride)
{
    
    // int idx = (blockIdx.x * blockDim.x + threadIdx.x + 1) * stride - 1;
    long long idx = (blockIdx.x * blockDim.x + threadIdx.x + 1) * stride - 1; // avoid integer overflow
    inputArray[idx] += inputArray[idx - (stride >> 1)];
}

__global__ void
downsweep_kernel(int N, int* inputArray, int stride)
{
    // int idx = ( blockIdx.x * blockDim.x + threadIdx.x + 1) * stride - 1;
    long long idx = ( blockIdx.x * blockDim.x + threadIdx.x + 1) * stride - 1; // avoid integer overflow
    int tmp = inputArray[idx - (stride >> 1)];
    inputArray[idx - (stride >> 1)] = inputArray[idx];
    inputArray[idx] += tmp;
}

然后是实现find_repeats功能,也不是很难,调用链是 fill_repeat_flags => exclusive_scan => scatter_flags,在实现第一个和第二个函数即可,没什么要注意的。

最后的得分如下,不清楚为什么数据量小的时候会有这么大差距

-------------------------
Scan Score Table:
-------------------------
-------------------------------------------------------------------------
| Element Count   | Ref Time        | Student Time    | Score           |
-------------------------------------------------------------------------
| 1000000         | 1.589           | 2.26            | 0.8788716814159293 |
| 10000000        | 14.072          | 12.157          | 1.25            |
| 20000000        | 27.321          | 22.288          | 1.25            |
| 40000000        | 53.277          | 42.531          | 1.25            |
-------------------------------------------------------------------------
|                                   | Total score:    | 4.628871681415929/5.0 |
-------------------------------------------------------------------------

-------------------------
Find_repeats Score Table:
-------------------------
-------------------------------------------------------------------------
| Element Count   | Ref Time        | Student Time    | Score           |
-------------------------------------------------------------------------
| 1000000         | 2.738           | 3.746           | 0.9136412172984516 |
| 10000000        | 20.158          | 20.511          | 1.25            |
| 20000000        | 39.049          | 38.149          | 1.25            |
| 40000000        | 76.606          | 70.088          | 1.25            |
-------------------------------------------------------------------------
|                                   | Total score:    | 4.663641217298451/5.0 |
-------------------------------------------------------------------------

render

实现一个简单的渲染器,refRenderer.cpp 文件已经给出了一个正确的串行方案,要求我们使用cuda语言实现并行版本。PA已经给出了一个可运行但是不正确的并行方案,它是以圆的颗粒度进行并行计算的,但是这破坏圆颜色之间的依赖关系,造成最终渲染结果的混乱。

__global__ void kernelRenderCircles() {

    int index = blockIdx.x * blockDim.x + threadIdx.x;

    if (index >= cuConstRendererParams.numCircles)  // index 是圆的索引
        return;
	// .....
    // for all pixels in the bonding box
    for (int pixelY=screenMinY; pixelY<screenMaxY; pixelY++) {
        float4* imgPtr = (float4*)(&cuConstRendererParams.imageData[4 * (pixelY * imageWidth + screenMinX)]);
        for (int pixelX=screenMinX; pixelX<screenMaxX; pixelX++) {
            float2 pixelCenterNorm = make_float2(invWidth * (static_cast<float>(pixelX) + 0.5f),
                                                 invHeight * (static_cast<float>(pixelY) + 0.5f));
            shadePixel(index, pixelCenterNorm, p, imgPtr);
            imgPtr++;
        }
    }
}
实现1

基于PA readme中提示,“There are two potential axes of parallelism in this assignment. One axis is parallelism across pixels another is parallelism across circles ”, 可以基于图像的每个像素进行并行计算,在kernel函数中再按照圆的依赖顺序依次遍历并调用shadePixel进行渲染。这样的实现是正确的,但是性能很差:

__global__ 
void kernelRenderPixels()
{
    int index = blockIdx.x * blockDim.x + threadIdx.x;
    int imageWidth = cuConstRendererParams.imageWidth;
    int imageHeight = cuConstRendererParams.imageHeight;

    if (index > imageWidth * imageHeight) 
    {
        return;
    }

    float invWidth = 1.f / imageWidth;
    float invHeight = 1.f / imageHeight;
    int pixelY = index / imageWidth;
    int pixelX = index % imageWidth;

    float2 pixelCenterNorm = make_float2(invWidth * (static_cast<float>(pixelX) + 0.5f),
                                         invHeight * (static_cast<float>(pixelY) + 0.5f));
    float4* imgPtr = (float4*)(&cuConstRendererParams.imageData[4 * index]);
    for (int circleIndex = 0; circleIndex < cuConstRendererParams.numCircles; ++circleIndex)
    {
        float3 circlePosition = *(float3*)(&cuConstRendererParams.position[3 * circleIndex]);
        shadePixel(circleIndex, pixelCenterNorm, circlePosition, imgPtr);
    }
}

CudaRenderer::render() {

    // 256 threads per block is a healthy number
    dim3 blockDim(256, 1);
    dim3 gridDim((image->width * image->height + blockDim.x - 1) / blockDim.x);
    kernelRenderPixels<<<gridDim, blockDim>>>();
    cudaDeviceSynchronize();
}

该实现的测试结果为:

--------------------------------------------------------------------------
| Scene Name      | Ref Time (T_ref) | Your Time (T)   | Score           |
--------------------------------------------------------------------------
| rgb             | 0.7622           | 0.6416          | 9               |
| rand10k         | 5.0429           | 77.2235         | 2               |
| rand100k        | 45.496           | 813.3614        | 2               |
| pattern         | 1.0921           | 8.9045          | 3               |
| snowsingle      | 29.2303          | 773.4755        | 2               |
| biglittle       | 27.4429          | 96.223          | 4               |
| rand1M          | 328.6162         | 8376.8588       | 2               |
| micro2M         | 606.1327         | 16774.7623      | 2               |
--------------------------------------------------------------------------
|                                    | Total score:    | 26/72           |
--------------------------------------------------------------------------
实现2

为了提升性能,我们把图片分成一个个16×16的小方块 ,每个threadblock负责这个小方块内的计算,主要包括:

  1. 并行地判断每个图片小方块是否与每个圆相交(并行判断256 个圆)
  2. 如果方块与圆相交,则再调用shadepixel进行渲染,如果不相交则直接跳过
#define BLOCKNUMX 16
#define BLOCKNUMY 16
#define BLOCKSIZE 256
__global__ 
void kernelRenderPixels()
{
    __shared__ int isBoxInCircle[BLOCKSIZE];

    int pixelX = blockIdx.x * blockDim.x + threadIdx.x;
    int pixelY = blockIdx.y * blockDim.y + threadIdx.y;

    int imageWidth = cuConstRendererParams.imageWidth;
    int imageHeight = cuConstRendererParams.imageHeight;

    float invWidth = 1.f / imageWidth;
    float invHeight = 1.f / imageHeight;

    if (pixelX >= imageWidth || pixelY >= imageHeight) 
    {
        return;
    }
    int boxL = blockIdx.x * blockDim.x;
    int boxR = (min(blockIdx.x * blockDim.x + blockDim.x, imageWidth));
    int boxB = blockIdx.y * blockDim.y;
    int boxT = (min(blockIdx.y * blockDim.y + blockDim.y, imageHeight));

    float boxLInv = boxL * invWidth;
    float boxRInv = boxR * invWidth;
    float boxTInv = boxT * invHeight;
    float boxBInv = boxB * invHeight;
    int linearThreadIndex =  threadIdx.y * blockDim.x + threadIdx.x;
    float2 pixelCenterNorm = make_float2(invWidth * (static_cast<float>(pixelX) + 0.5f),
                                        invHeight * (static_cast<float>(pixelY) + 0.5f));
    for (int batchStartIndexForCircles = 0; 
         batchStartIndexForCircles < cuConstRendererParams.numCircles;
         batchStartIndexForCircles += BLOCKSIZE)
    {
        int indexForCircles = batchStartIndexForCircles + linearThreadIndex;
        if (indexForCircles >= cuConstRendererParams.numCircles)
        {
            isBoxInCircle[linearThreadIndex] = 0;
        }
        else 
        {
            float circleX = cuConstRendererParams.position[3 * indexForCircles];
            float circley = cuConstRendererParams.position[3 * indexForCircles + 1];
            float circleRadius = cuConstRendererParams.radius[indexForCircles];
            isBoxInCircle[linearThreadIndex] = circleInBoxConservative(circleX, circley, 
                                                           circleRadius, boxLInv, boxRInv, boxTInv, boxBInv) ? circleInBox(circleX, circley, 
                                                           circleRadius, boxLInv, boxRInv, boxTInv, boxBInv) : 0;
        }
        __syncthreads();


        for (int i = batchStartIndexForCircles; i< batchStartIndexForCircles + BLOCKSIZE; ++i)
        {
            if (isBoxInCircle[i % BLOCKSIZE])
            {
                float4* imgPtr = (float4*)(&cuConstRendererParams.imageData[4 * (pixelY * imageWidth + pixelX)]);
                float3 circlePosition = *(float3*)(&cuConstRendererParams.position[3 * i]);
                shadePixel(i, pixelCenterNorm, circlePosition, imgPtr);
            }
        }
        __syncthreads();
    }

}

void
CudaRenderer::render() 
{
    dim3 blockDim(BLOCKNUMX, BLOCKNUMY);
    dim3 gridDim((image->width + blockDim.x - 1) / blockDim.x, (image->height + blockDim.y - 1) / blockDim.y);
    kernelRenderPixels<<<gridDim, blockDim>>>();  
}

注意,最后一个__syncthreads是必须得加的,否则会导致线程安全问题(为了保护某一cuda线程在一次for循环中,isBoxInCircle不被复写)。

最后再测下性能:

--------------------------------------------------------------------------
| Scene Name      | Ref Time (T_ref) | Your Time (T)   | Score           |
--------------------------------------------------------------------------
| rgb             | 0.8281           | 1.5998          | 6               |
| rand10k         | 5.1915           | 37.4328         | 3               |
| rand100k        | 45.6692          | 345.5081        | 3               |
| pattern         | 1.0687           | 4.5626          | 4               |
| snowsingle      | 28.8797          | 295.1067        | 2               |
| biglittle       | 26.8013          | 60.9962         | 6               |
| rand1M          | 339.1135         | 3399.8742       | 2               |
| micro2M         | 609.6217         | 6839.2456       | 2               |
--------------------------------------------------------------------------
|                                    | Total score:    | 28/72           |
--------------------------------------------------------------------------

Amazing!可以说是没有任何的性能提升呢......

实现3

没法了,去参考了别人的实现,发现他们除了对图像进行分块,还在第一个循环内使用sharedMemExclusiveScan进行优化:

__global__ 
void kernelRenderPixels()
{
    __shared__ uint isBoxInCircle[BLOCKSIZE];
    __shared__ uint prefixSumOutput[BLOCKSIZE];
    __shared__ uint prefixSumScratch[2 * BLOCKSIZE];
    __shared__ int inBoxCircleIndexes[BLOCKSIZE];
	// ...
   
    int linearThreadIndex =  threadIdx.y * blockDim.x + threadIdx.x;
    float2 pixelCenterNorm = make_float2(cuConstRendererParams.invWidth * (static_cast<float>(pixelX) + 0.5f),
                                        cuConstRendererParams.invHeight * (static_cast<float>(pixelY) + 0.5f));
    isBoxInCircle[linearThreadIndex] = 0;
    inBoxCircleIndexes[linearThreadIndex] = -1;
    for (int batchStartIndexForCircles = 0; 
         batchStartIndexForCircles < cuConstRendererParams.numCircles;
         batchStartIndexForCircles += BLOCKSIZE)
    {
        int indexForCircles = batchStartIndexForCircles + linearThreadIndex;
        if (indexForCircles < cuConstRendererParams.numCircles)
        {
            float circleX = cuConstRendererParams.position[3 * indexForCircles];
            float circley = cuConstRendererParams.position[3 * indexForCircles + 1];
            float circleRadius = cuConstRendererParams.radius[indexForCircles];
            isBoxInCircle[linearThreadIndex] =  circleInBox(circleX, circley, circleRadius, boxLInv, boxRInv, boxTInv, boxBInv);
        }
        __syncthreads();

        // shoudl use sharedMemExclusiveScan to improve permance
        // but why the performance is ** far far** better than the one that do not do the execlusive scan improvement?
        sharedMemExclusiveScan(linearThreadIndex, isBoxInCircle, prefixSumOutput, prefixSumScratch, BLOCKSIZE);

        if (isBoxInCircle[linearThreadIndex]) {
            inBoxCircleIndexes[prefixSumOutput[linearThreadIndex]] = indexForCircles;
        }

        __syncthreads();

        int numOfIntescetedCircles = prefixSumOutput[BLOCKSIZE - 1] + isBoxInCircle[BLOCKSIZE - 1];
        for (int i = 0; i < numOfIntescetedCircles; ++i) // 只循环遍历实际相交的圆
        {
            float4* imgPtr = (float4*)(&cuConstRendererParams.imageData[4 * (pixelY * imageWidth + pixelX)]);
            float3 circlePosition = *(float3*)(&cuConstRendererParams.position[3 * inBoxCircleIndexes[i]]);
            shadePixel(inBoxCircleIndexes[i], pixelCenterNorm, circlePosition, imgPtr);
        }
    }

起初我以为这个优化是微不足道的,因为第二层for循环最多也就循环256次(我认为,这和和测试集圆的数量——100k比起来可以约等于O(1)的复杂度),调用sharedMemExclusiveScan带来的cost可能还大于减少循环的收益。但测试结果显示,非常的Amazing啊,优化了非常多:

--------------------------------------------------------------------------
| Scene Name      | Ref Time (T_ref) | Your Time (T)   | Score           |
--------------------------------------------------------------------------
| rgb             | 0.6982           | 0.702           | 9               |
| rand10k         | 6.4562           | 7.4349          | 9               |
| rand100k        | 45.3627          | 55.6973         | 8               |
| pattern         | 1.0167           | 0.9699          | 9               |
| snowsingle      | 34.3391          | 37.9076         | 9               |
| biglittle       | 25.9874          | 53.0712         | 6               |
| rand1M          | 322.1914         | 325.6942        | 9               |
| micro2M         | 586.4424         | 597.0768        | 9               |
--------------------------------------------------------------------------
|                                    | Total score:    | 68/72           |
--------------------------------------------------------------------------

但是为什么呢?我只能猜测,实现3的for循环内部没有了if循环的分支,大大减少了execution divergence

// 实现2
for (int i = batchStartIndexForCircles; i< batchStartIndexForCircles + BLOCKSIZE; ++i)
{
    if (isBoxInCircle[i % BLOCKSIZE])
    {
        float4* imgPtr = (float4*)(&cuConstRendererParams.imageData[4 * (pixelY * imageWidth + pixelX)]);
        float3 circlePosition = *(float3*)(&cuConstRendererParams.position[3 * i]);
        shadePixel(i, pixelCenterNorm, circlePosition, imgPtr);
    }
}
// 实现3
int numOfIntescetedCircles = prefixSumOutput[BLOCKSIZE - 1] + isBoxInCircle[BLOCKSIZE - 1];
for (int i = 0; i < numOfIntescetedCircles; ++i) // 只循环遍历实际相交的圆
{
    float4* imgPtr = (float4*)(&cuConstRendererParams.imageData[4 * (pixelY * imageWidth + pixelX)]);
    float3 circlePosition = *(float3*)(&cuConstRendererParams.position[3 * inBoxCircleIndexes[i]]);
    shadePixel(inBoxCircleIndexes[i], pixelCenterNorm, circlePosition, imgPtr);
}

本想用ncu等性能测试工具验证这个猜想的,但是捣鼓了一晚上还是没对齐ncu、diver和cudaruntime的版本,所以就先这样吧。

此外,先调circleInBoxConservative再调circleInBox的方法,没有做到有效优化;把第一层for循环的if去掉了,也没什么优化:

    for (int batchStartIndexForCircles = 0; 
         batchStartIndexForCircles < cuConstRendererParams.numCircles;
         batchStartIndexForCircles += BLOCKSIZE)
    {
        int indexForCircles = batchStartIndexForCircles + linearThreadIndex;

        float circleX = cuConstRendererParams.position[3 * indexForCircles];
        float circley = cuConstRendererParams.position[3 * indexForCircles + 1];
        float circleRadius = cuConstRendererParams.radius[indexForCircles];
        isBoxInCircle[linearThreadIndex] =  circleInBox(circleX, circley, circleRadius, boxLInv, boxRInv, boxTInv, boxBInv);
        __syncthreads();
        // ...

Assignment 4

环境 配置使用conda 安装对应pytorch版本 [2.1.2]( https://pytorch.org/get-started/previous-versions/#v212:~:text=org/whl/cpu- ,v2.1.2,-Conda)

但还需解决一些版本问题和依赖安装问题:

  • 改mkl版本: conda install mkl=2023.1.0
  • 改numpy版本 :conda install numpy=1.26 -y
  • 安装低版本的 setuptools: pip install setuptools==59.6.0 --force-reinstall
  • 安装依赖 tiktoken: conda install tiktoken
  • 用 conda list > env.txt 命令导出了依赖,可以参考下

整个实验就是对矩阵乘法的优化,至于GPT、大模型原理啥的另请高明吧,但做这个实验也不太需要知道这些原理。最后一个实验没学到啥的感觉,part4也跳过了。

Part1 A Simple (But Not So Efficient) Implementation of Attention

按照Readme中给出的步骤依次实现即可

    for (int b = 0; b < B; ++b)
    {
        for (int h = 0; h < H; ++h)
        {
            for (int i = 0; i < N; ++i)
            {
                for (int j = 0; j < N; ++j)
                {
                    float sum = 0.0f;
                    for (int k = 0; k < d; ++k)
                    {
                        sum += fourDimRead(Q, b, h, i, k, H, N, d) * fourDimRead(K, b, h, j, k, H, N, d);
                    }
                    twoDimWrite(QK_t, i, j, N, sum);
                }
            }

            for (int i = 0; i < N; ++i) 
            {
                float sum = 0.0;
                for (int j = 0; j < N; ++j) 
                {
                    sum += std::exp(twoDimRead(QK_t, i, j, N));
                }

                for (int j = 0; j < N; ++j) 
                {
                    float val = std::exp(twoDimRead(QK_t, i, j, N)) / sum;
                    twoDimWrite(QK_t, i, j, N, val);
                }
            }

            for (int i = 0; i < N; i++) 
            {
                for (int j = 0; j < d; j++) 
                {
                    float sum = 0.0;
                    for (int k = 0; k < N; k++) 
                    {
                        sum += twoDimRead(QK_t, i, k, N) * fourDimRead(V, b, h, k, j, H, N, d);
                    }
                    fourDimWrite(O, b, h, i, j, H, N, d, sum);
                }
            }
        }
    }

Part2 Blocked Matrix Multiply and Unfused Softmax

基于block进行优化,根据课程ppt上的该就行了。理论上BlockSize的大小,要和cacheline一致才能达到最大的加速比。查看cahceline大小:

cat /sys/devices/system/cpu/cpu1/cache/index0/coherency_line_size 64

所以把BLOCK_SIZE 定为16了

实现如下:

  for (int b = 0; b < B; ++b)
    {
        for (int h = 0; h < H; ++h)
        {
            for (int i_block = 0; i_block < N; i_block += BLOCK_SIZE)
            {
                for (int j_block = 0; j_block < N; j_block += BLOCK_SIZE)
                {
                    for (int k_block = 0; k_block < d; k_block += BLOCK_SIZE)
                    {

                        for (int i = 0; i < BLOCK_SIZE; ++i)
                        {
                            if (i + i_block >= N) break;
                            for (int j = 0; j < BLOCK_SIZE; ++j)
                            {
                                if (j + j_block >= N) break;
                                float block_sum = 0.f;
                                for (int k = 0; k < BLOCK_SIZE; ++k)
                                {
                                    if (k + k_block >= d) break;
                                    block_sum += fourDimRead(Q, b, h, i + i_block, k + k_block, H, N, d) * fourDimRead(K, b, h, j + j_block, k + k_block, H, N, d);
                                }
                                float origin = twoDimRead(QK_t, i + i_block, j + j_block, N);
                                twoDimWrite(QK_t, i + i_block, j + j_block, N, block_sum + origin);
                            }
                        }
                    }
                }
            }

            for (int i = 0; i < N; ++i) 
            {
                float sum = 0.0;
                for (int j = 0; j < N; ++j) 
                {
                    sum += std::exp(twoDimRead(QK_t, i, j, N));
                }

                for (int j = 0; j < N; ++j) 
                {
                    float val = std::exp(twoDimRead(QK_t, i, j, N)) / sum;
                    twoDimWrite(QK_t, i, j, N, val);
                }
            }

            for (int i_block = 0; i_block < N; i_block += BLOCK_SIZE)
            {
                for (int j_block = 0; j_block < d; j_block += BLOCK_SIZE)
                {
                    for (int k_block = 0; k_block < N; k_block += BLOCK_SIZE)
                    {

                        for (int i = 0; i < BLOCK_SIZE; ++i)
                        {
                            if (i + i_block >= N) break;
                            for (int j = 0; j < BLOCK_SIZE; ++j)
                            {
                                if (j + j_block >= d) break;
                                float block_sum = 0.f;
                                for (int k = 0; k < BLOCK_SIZE; ++k)
                                {
                                    if (k + k_block >= N) break;
                                    block_sum += twoDimRead(QK_t, i + i_block, k + k_block, N) * fourDimRead(V, b, h, k + k_block, j + j_block, H, N, d);
                                }
                                float origin = fourDimRead(O, b, h, i + i_block, j + j_block, H, N, d);
                                fourDimWrite(O, b, h, i + i_block, j + j_block, H, N, d, block_sum + origin);
                            }
                        }
                    }
                }
            }

        }
    }

结果如下所示,有时候程序性能比ref程序好,但有时候也会差很多。

-----RUNNING REFERENCE IMPLEMENTATION-----
STUDENT - BLOCKED MATMUL + UNFUSED SOFTMAX statistics
cpu time:  196.303ms
mem usage:  4718592 bytes
-----RUNNING STUDENT IMPLEMENTATION-----
STUDENT - BLOCKED MATMUL + UNFUSED SOFTMAX statistics
cpu time:  185.333ms
mem usage:  4718592 bytes

参考了其他人的实现后发现,把if分支去掉后能得到大幅提升:

    for (int b = 0; b < B; ++b)
    {
        for (int h = 0; h < H; ++h)
        {
            for (int i_block = 0; i_block < N; i_block += BLOCK_SIZE)
            {
                for (int j_block = 0; j_block < N; j_block += BLOCK_SIZE)
                {
                    for (int k_block = 0; k_block < d; k_block += BLOCK_SIZE)
                    {
                        int i_limit = std::min(N - i_block, BLOCK_SIZE);
                        int j_limit = std::min(N - j_block, BLOCK_SIZE);
                        int k_limit = std::min(d - k_block, BLOCK_SIZE);
                        for (int i = 0; i < i_limit; ++i)
                        {
                            for (int j = 0; j < j_limit; ++j)
                            {
                                float block_sum = 0.f;
                                for (int k = 0; k < k_limit; ++k)
                                {
                                    block_sum += fourDimRead(Q, b, h, i + i_block, k + k_block, H, N, d) * fourDimRead(K, b, h, j + j_block, k + k_block, H, N, d);
                                }
                                float origin = twoDimRead(QK_t, i + i_block, j + j_block, N);
                                twoDimWrite(QK_t, i + i_block, j + j_block, N, block_sum + origin);
                            }
                        }
                    }
                }
            }

            for (int i = 0; i < N; ++i) 
            {
                float sum = 0.0;
                for (int j = 0; j < N; ++j) 
                {
                    sum += std::exp(twoDimRead(QK_t, i, j, N));
                }

                for (int j = 0; j < N; ++j) 
                {
                    float val = std::exp(twoDimRead(QK_t, i, j, N)) / sum;
                    twoDimWrite(QK_t, i, j, N, val);
                }
            }

            for (int i_block = 0; i_block < N; i_block += BLOCK_SIZE)
            {
                for (int j_block = 0; j_block < d; j_block += BLOCK_SIZE)
                {
                    for (int k_block = 0; k_block < N; k_block += BLOCK_SIZE)
                    {

                        int i_limit = std::min(N - i_block, BLOCK_SIZE);
                        int j_limit = std::min(d - j_block, BLOCK_SIZE);
                        int k_limit = std::min(N - k_block, BLOCK_SIZE);
                        for (int i = 0; i < i_limit; ++i)
                        {
                            for (int j = 0; j < j_limit; ++j)
                            {
                                float block_sum = 0.f;
                                for (int k = 0; k < k_limit; ++k)
                                {
                                    block_sum += twoDimRead(QK_t, i + i_block, k + k_block, N) * fourDimRead(V, b, h, k + k_block, j + j_block, H, N, d);
                                }
                                float origin = fourDimRead(O, b, h, i + i_block, j + j_block, H, N, d);
                                fourDimWrite(O, b, h, i + i_block, j + j_block, H, N, d, block_sum  + origin);
                            }
                        }
                    }
                }
            }

        }
    }


-----RUNNING REFERENCE IMPLEMENTATION-----
REFERENCE - BLOCKED MATMUL + UNFUSED SOFTMAX statistics
cpu time:  164.46ms
mem usage:  4718592 bytes
-----RUNNING STUDENT IMPLEMENTATION-----
STUDENT - BLOCKED MATMUL + UNFUSED SOFTMAX statistics
cpu time:  131.369ms
mem usage:  4718592 bytes

其中的原因可能有两个:

  • CPU对if分支有预测,如果预测错误整条流水线清空,浪费大量CPU时间
  • if分支越多,编译器越难对代码进行SIMD处理,导致性能下降

4个实验下来,这里是第二次看到消除if分支带来的性能提升(第一次有较大体会的是asst3)

Part3 Fused Attention

回顾part2,可以看到softmax不能进行Block优化,但可以被fuse,也就是说Q * K^t 的一行算完后,就进行softmax操作,然后计算PV的一行并存储到结果中的一行中。这样做我们仅需要一维的中间变量,对内存使用量有一定的优化,而且这样的算法中每一个batch、每一个head、每一行的for循环计算是互不依赖的,这意味着我们可以用多线程使其并行化。试验要求我们使用OPENMP进行优化。代码如下:

    #pragma omp parallel for collapse(3)
    for (int b = 0; b < B; ++b)
    {
        for (int h = 0; h < H; ++h)
        {
            for (int i = 0; i < N; ++i)
            {
                // YRow is moved inside so each OpenMP thread gets a local copy.
                at::Tensor ORowTensor = temp.index({torch::indexing::Slice(
                    omp_get_thread_num(), torch::indexing::None)});
                std::vector<float> ORow = formatTensor(ORowTensor);

                for (int j = 0; j < N; ++j)
                {
                    float sum = 0.0f;
                    for (int k = 0; k < d; ++k)
                    {
                        sum += fourDimRead(Q, b, h, i, k, H, N, d) * fourDimRead(K, b, h, j, k, H, N, d);
                    }
                    ORow[j] = sum;
                }

                float sum = 0.0;
                for (int j = 0; j < N; ++j) 
                {
                    ORow[j] = std::exp(ORow[j]);
                    sum += ORow[j];
                }
                for (int j = 0; j < N; ++j) 
                {
                    ORow[j] /= sum;
                }

                for (int j = 0; j < d; ++j) 
                {
                    float sum = 0.0;
                    for (int k = 0; k < N; ++k) 
                    {
                        sum += ORow[k] * fourDimRead(V, b, h, k, j, H, N, d);
                    }
                    fourDimWrite(O, b, h, i, j, H, N, d, sum);
                }
            }
        }
    }

测试结果如下:

-----RUNNING REFERENCE IMPLEMENTATION-----
REFERENCE - FUSED ATTENTION statistics
cpu time:  38.951ms
mem usage:  557056 bytes
-----RUNNING STUDENT IMPLEMENTATION-----
Self CPU time total: 35.387ms
STUDENT - FUSED ATTENTION statistics
cpu time:  35.336ms
mem usage:  557056 bytes
posted @ 2026-06-28 22:46  别杀那头猪  阅读(4)  评论(0)    收藏  举报