《并行计算与高性能计算》简记

并行计算与高性能计算》阅读简记。

第一部分并行计算介绍

并行计算的基本定律：
- Amdahl定律
- 强标度：总体问题规模确定，处理器数量增加，求解时间的变化。
- 弱标度：单个处理器处理的问题规模确定，处理器数量增加，求解时间的变化。
并行方法(模型)：共享内存，消息传输，数据并行，流处理(GPU)。
算数强度：每个内存(4/8字节的数据)执行浮点运算次数。
理论性能峰值：虚拟核心数 × 频率 × 向量位 × FMA / (32或64)。
理论内存带宽：传输速率 × 通道数 × 访问字节数 × 插槽数。
实证测试：stream、cpufp。
Roofline模型

面向内存带宽的数据结构设计：
- 连续的二维数组分配。
- 结构数组、数组结构、混合结构，根据数据的使用情况确定合适的方式。
除了尽量将连续访问的内存分配在一起，一次完成可预知大小的分配比碎片化的分配要效果更好

第二部分 CPU并行

向量化的部分介绍了手动向量化和编译器自动向量化的方式，比较常规。

OpenMP多线程并行

OpenMP的部分除了介绍基本用法外，还介绍了一些重要的影响性能的细节：

first touch：内存位置会分配到在第一次访问该地址的线程所在的位置，在此之前仅存在于虚拟地址空间的页表中，因此让每个线程first touch需要计算的数据有利于访存
高级OpenMP：当程序有多次OpenMP并行时，可以只启动并行段一次，将整体程序全部包含在其中，串行的部分通过线程号控制主线程处理，并且通过nowait最小化barrier和可能的同步操作

书中给出了几个示例。这里也尝试一下前面提到的内存相关的优化效果。

Example: SplitStencil

opt	time
baseline	27.459721
memopt	14.876064
omp4线程	10.687703
omp(一次并行启动4线程)	10.937465

这里一次启动没有什么效果，暂不确定原因。但是memopt作用很明显。另外多线程效果也非常差。在另一台机器上多线程效果好一些，一次启动同样没有明显的效果，可能是启动开销并不大，memopt也没有效果，应该是原本的分配方式也刚好malloc到了统一的内存位置。关于性能上的差异和产生的问题，需要进一步profile才能确定原因。

Example: PrefixScan

排他前缀和的并行算法，首先将数据整体进行分块，每个线程首先计算自己块内的排他前缀和，此时每个线程的第一个位置为0，然后由主线程串行更新偏移，设定每个线程j的第一个位置为前j-1个线程的排他前缀和，这个值为j-1的起始位置+线程j-1的最后一个位置的排他前缀和+最后一个位置的值，特殊处理线程块内只有两个元素的情况，此时上一区间没有排他前缀和或者说为0，已经在处理上一个线程数据时修改了。最后再并行将每个块的第一个值也就是真正的前缀加到块内每一个数据上。

void PrefixScan (int *input, int *output, int length) 
{
   // Get the total number of threads and thread_id
   int nthreads = 1;
   int thread_id   = 0;
#ifdef _OPENMP
   nthreads = omp_get_num_threads();
   thread_id = omp_get_thread_num();
#endif

   // Compute the range for which this thread is responsible.
   int tbegin = length * ( thread_id     ) / nthreads;
   int tend   = length * ( thread_id + 1 ) / nthreads;

   // Only perform this operation if there is a positive number of entries.
   if ( tbegin < tend ) {
       // Do an exclusive scan for each thread
       output[tbegin] = 0;
       for ( int i = tbegin + 1 ; i < tend ; i++ ) {
          output[i] = output[i-1] + input[i-1];
       }
   }
   if (nthreads == 1) return;

   // Do adjustment to prefix scan for the beginning value for each thread
#ifdef _OPENMP
   // Wait until all threads get here. 
#pragma omp barrier
    
   // On the master thread compute the beginning offset for each thread
   if (thread_id == 0) {
      for ( int i = 1 ; i < nthreads ; i ++ ) {
         int ibegin = length * ( i - 1 ) / nthreads;
         int iend   = length * ( i     ) / nthreads;
         //ibegin是前一个线程块之前的数据的前缀和，iend是当前块的第一个位置
         if ( ibegin < iend ) 
            output[iend] = output[ibegin] + input[iend-1];

         if ( ibegin < iend - 1 )
            output[iend] += output[iend-1];
      }
   }
#pragma omp barrier
   // Start all threads again
   // Apply the offset to the range for this thread.
#pragma omp simd
   for ( int i = tbegin + 1 ; i < tend ; i++ ) {
      output[i] += output[tbegin];
   }
#endif
}